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Preface 


In  recenc  years , /advances  in  composite  materials  and  active  control 
principles  have  lead  to  renewed  interest  in  the  Forward  Swept  Wing^as  a 
feasible  aircraft  design  alternative.  Sweeping  a  wing  forward  results  in 
low  static  divergence  speeds  and  additional  aeroelastic  instabilities 

which  must  be  adequately  controlled  foraa  aircraft  employing  this  wing 

ttj  %!  i  £  tS  ff  S  to  U'Vt  isyr  t  s  J 


l>  .•*  ^  s  k  s  to  %-fr  a 

design  to  have  an  acceptable  flight  envelope.  ^Contained  herein  is  a  study 
aimed  at  dcmouBtrating^the  utility  of  applying  active  feedback  control 
principles  to  suppress  ther  aeroelastic  instabilities  associated  with  a  h  F5U/ 
Forward  Swept  Wing  ‘desirn.  -s  Analytical  techniques  are  presented  that  are 

useful  in  the  analysis  and  Synthesis  of  active  control  laws  using  optimal 

Ss — v>  <3/'V  ^  /  1  J*  I X  ) 

control  theory  methodolgy.  '  • 0  J 
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Abstract 


/Analytical  studies  Were  conducted"^  investigate) the  potential  of 

'  r-  ^  ■ 

l  / 

applying  optimal  control  theory  techniques  to  synthesis  of  active 
flutter  suppression  control  lavs.  For  an  example  application,  arForwardV 


fuselage  model  previously  analyzed  by  Thomas  Noll  of  -the  Air 


^-^^oree- 


,  U  Z,j!  ’  O— 

Flight  Dynamics  Laboratory  was  utilized;  Through ‘'the  use  of  Pade' 


approximants  to  represent 


unsteady 


aerodynamic  forces,  the  equations 


of  motion  are  written  in  standard  state  space  form.  Linear  optimal  reg¬ 
ulator  theory  is  then  applied  to  determine  particular  sets  of  gains  which 

minimize  a  quadratic  cost  function  in  terms  of  the  states  and  controls. 

r> 

The  control  laws  are  developed  at  a  design  flight  condition  which  in- 

•»> 

creases  the  onset  of  the  lowest  instability  speed  20Z  above  the  wing 
bending/torsion  Instability  speed.  The  optimal  control  law  is  then  ap¬ 


plied  at  off-design  flight  conditions  to  assess  the  robustness  of  the  opt¬ 
imal  regulator.  r.  \  ,v  «’■  X 
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ACTIVE  SUPPRESSION  OF  AEROELASTIC  INSTABILITIES 
ON  A  FORWARD  SWEPT  WING  OSIN G  A 
LINEAR  OPTIMAL  REGULATOR 


I.  Introduction 


In  recent  years,  the  forward  swept  wing  has  gained  renewed  Interest 
among  the  aerospace  community  as  being  a  feasible  aircraft  design  altern¬ 
ative.  This  is  not  a  new  concept  as  the  forward  swept  wing  has  long  been 
recognized  as  being  able  to  provide  improved  performance  benefits  over 
conventional  aft  swept  wings,  provided  potential  aeroelastic  problems 
could  be  solved.  The  aeroelastic  instabilities  associated  with  the  oper¬ 
ation  of  a  forward  swept  wing  Include  instabilities  not  typically  encount¬ 
ered  within  the  operational  flight  envelope  of  conventional  aft  swept  wing 
aircraft,  such  as  static  divergence  and  body  freedom  flutter.  Convention¬ 
al  passive  methods  for  preventing  these  instabilities  typically  result  in 
significant  performance  penalties  to  the  aircraft,  negating  the  Improved 
performance  gained  from  using  a  forward  swept  wing  design.  Therefore, 
there  is  considerable  Interest  in  developing  other  methods  of  preventing 
aeroelastic  instabilities  (or  increasing  the  minimum  airspeed  at  which 
they  occur)  that  can  be  used  in  place  of,  or  in  combination  with,  passive 
methods. 

A  method  being  considered  for  preventing  aeroelastic  instabilities 
Involves  the  utilization  of  an  active  feedback  control  system.  An  act- 
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ive  control  system  utilizes  an  aerodynamic  control  surface,  or  surfaces, 
commanded  by  signals  through  an  an  appropriate  control  law.  The  funda¬ 
mental  principles  behind  the  use  of  this  technique  have  been  well  doc¬ 
umented  as  a  result  of  the  significant  amount  of  research  conducted  in 
the  1960s  and  1970s  to  develop  active  flutter  suppression  technology 
(Ref  1-4),  and  in  the  early  1980s  to  advance  adaptive  control  principles. 
With  these  advances  in  active  control  technology,  an  active  flutter  sup¬ 
pression  control  system  is  now  a  feasible  solution  to  the  aeroelastlc 
problems  associated  with  forward  swept  wing  aircraft,  while  offering 
significantly  less  weight  and  performance  penalties  than  alternative 
solutions. 

There  have  been  several  studies  conducted  in  recent  years  (Ref  5-6) 
pertaining  to  the  use  of  classical  control  theory  methods  for  active 
flutter  suppression  on  aircraft  employing  a  forward  swept  wing  design. 
Griffen  and  Eastep  initially  investigated  the  use  of  a  simple  active 
feedback  control  system  to  control  divergence  and  bending/torsion  flut¬ 
ter  separately  on  a  cantilever  wing  configuration.  At  the  time  of  this 
study,  static  divergence  was  believed  to  be  the  most  critical  aeroelastlc 
Instability  but  when  a  wing/fuselage  model  was  tested  free  in  pitch 
(Ref  7-8),  body  freedom  flutter  was  discovered  to  occur  at  a  velocity 
lower  than  the  static  divergence  speed.  This  instability  is  a  result  of 
coupling  between  the  airplane  rigid  body  pitch  mode  (short  period)  and 
the  wing  first  bending  mode.  A  more  recent  study  conducted  by  Noll  et  al 
(Ref  9)  was  directed  at  developing  control  laws  for  a  cantilever  flexible 
wing  and  for  a  wing/fuselage  model  free  in  pitch.  In  this  study,  Noll 
developed  a  multiple  input-multiple  output  (MIMO)  control  law  by  treat- 


log  Che  system  as  two  single  input-single  output  (SISO)  sub-systens,  and 
applying  classical  control  theory  methodology.  Two  active  control  sur¬ 
faces  were  utilized  along  with  appropriate  feedback  compensation,  thereby 
making  the  control  law  synthesis  an  iterative  procedure.  While  his  final 
control  law  design  did  accomplish  the  objectives  of  the  study,  it  will  be 
demonstrated  herein  that  rather  straight  forward,  optimal  control  theory 
techniques  can  be  utilized  for  the  control  lav  synthesis,  thereby  elim¬ 
inating  the  iterative  design  process.  Since  in  actuality  most  aircraft 
control  problems  are  MIMO  systems,  optimal  control  theory  techniques 
would  seem  rather  well  suited  for  active  control  law  synthesis. 

To  apply  optimal  regulator  theory  for  the  synthesis  of  an  active 
flutter  suppression  control  law,  it  is  first  necessary  to  obtain  the 
aeroelastlc  equations  of  motion  in  the  form  of  constant  coefficient  dif¬ 
ferential  equations.  The  problem  arises  in  representing  the  unsteady 
aerodynamics  for  arbitrary  motion  in  constant  coefficient  differential 
equation  form  so  as  to  be  incorporated  into  the  equations  of  notion 
easily.  In  the  study  conducted  by  Noll,  the  generalized  aerodynamic 
forces  for  oscllllatory  motion  were  computed  using  the  doublet-lattice 
method  of  Ref  10.  To  compute  the  unsteady  aerodynamics  for  arbitrary 
motion,  each  generalized  force  was  represented  by  a  Pade'  approxlmant , 
that  is,  a  ratio  of  polynomials  in  the  Laplace  variable  s.  This  enabled 
the  unsteady  aerodynamic  forces  to  be  easily  incorporated  into  the  trans¬ 
formed  equations  of  motion  to  be  used  as  a  foundation  for  the  control 
law  synthesis  procedure. 


The  purpose  of  this  thesis  is  to  present  some  analytical  techniques 
that  are  useful  for  the  analysis  and  synthesis  of  active  flutter  suppres- 


sion  control  laws  using  optimal  regulator  theory.  Control  laws  were 
synthesized  for  Noll's  cantilever  forward  swept  wing  model  at  a  design 
flight  condition  which  increases  the  onset  of  the  lowest  instability 
speed  20%  above  the  wing  bending/torsion  instability  speed.  The  optimal 
control  law  is  then  applied  at  off-design  flight  conditions  to  assess  the 
robustness  of  the  optimal  regulator.  The  effect  of  varying  the  state 
weighting  matrix  at  the  design  flight  condition  is  also  assessed  along 
with  the  adequacy  of  developing  the  control  law  based  on  a  reduced  order 
system  model  which  neglects  the  high  frequency  aerodynamic  lag  states. 


Gl 
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II.  Theoretical  Developement 


a 


Selected  Model  Configuration 

It  is  envisioned  that  a  fighter-type  aircraft  employing  a  forward 
swept  wing  design  could  evolve  to  the  point  of  operational  deployment, 
whereas  it  would  be  expected  to  exhibit  operational  characteristics  sim¬ 
ilar  to  those  of  conventional  aft  swept  wing  fighters.  Typically,  a 
fighter  aircraft  exhibiting  a  multi-role  capability  will  carry  externally 
mounted  stores  both  on  the  fuselage  and  uoder  the  wings.  However,  the 
adverse  mass  and  lnertls  distribution  on  the  wings  caused  by  the  external 
stores  have  traditionally  resulted  in  bending/torsion  flutter  problems 
occurring  within  the  operational  flight  envelope.  For  these  reasons,  the 
model  configuration  chosen  by  Noll  was  characterized  as  having  a  clamped 
divergence  instability  in  close  proximity  to  a  classical  bendlng/torslon 
flutter  mode,  typical  of  what  would  occur  for  a  critical  external  store 
configuration.  Figure  1  contains  a  schematic  of  the  model  selected  for 
the  current  analyses  and  shows  the  key  details  and  dimensions,  and  the 
relative  sizes  of  the  components. 

To  calculate  the  natural  frequencies  and  mode  shapes  of  the  forward 
swept  wing  configuration,  Noll  developed  a  finite  element  representation 
of  the  model  for  use  with  the  NASTRAN  program.  For  the  cantilever  anal¬ 
yses,  the  bar  representing  the  fuselage  was  restrained  in  such  a  manner 
that  no  motion  was  permitted  along  the  bar  or  wing  root.  The  elastic 
modes  used  in  the  analyses  of  this  study  included  the  first  two  bending 
modes  and  the  first  torsion  mode  of  the  wing.  It  was  determined  that 
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Figure  1.  Planform  of  Forward  Swept  Wing  Model 


higher  frequency  nodes  could  be  eliulnated  when  flutter  analyses  showed 
that  they  had  no  appreciable  affect  on  the  instabilities  of  Interest. 


State  Space  Equations  of  Motion 

To  apply  optimal  control  theory  techniques  for  control  law  synthesis. 
It  is  desired  to  obtain  the  aeroelastic  equations  of  notion  as  a  set  of 
linear  first  order  differential  equations  of  the  fora 

X(t)  -  (A]X(t)  +  (B)u(t)  (2-1) 
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where  X(t)  it  a  real  n-dieenaional  column  vector  containing  the  “states" 
of  the  system,  and  u(t)  Is  a  real  m-dimensional  column  vector  which  is  the 
input  control  vector.  The  equations  of  motion  of  a  flexible  aircraft  with 
control  surfaces  present  can  be  represented  as 

[M]£<t)  +  [C]£(t)  +  [K]£(t)  +  l/2,V2s[Q(k)]£(t) 

+  lMc]£.(t)  +  l/2pV2?[Qc (101^(0  -  0  (2-2) 

where  the  generalised  unsteady  aerodynamic  forces  are  contained  in  the<  ■ 
[Q(k) ]  and  [Qc(k)]  matrices  for  the  wing  and  control  surfaces  respect¬ 
ively.  A  more  detailed  derivation  of  the  equations  of  motion  and  an  ex¬ 
planation  of  the  aerodynamic  modeling  used  in  the  current  analysis  is 
contained  in  Appendix  A  and  further  detail  can  also  be  obtained  from  Ref 
9.  When  the  Laplace  Transform  of  Eq  2-2  is  taken,  substitutions  for  the 
unsteady  aerodynamics  are  made  and  powers  of  s  are  equated  with  time  der¬ 
ivatives,  resulting  in  the  equations  of  motion  being  represented  as: 

£(t)  +  [a  j£(t)  +  la^jtft)  +  la^sJt)  +  [a  J^(t)  - 

ib3J£(t)  +  lb2]i*cU)  +  lV*c(t)  +  lV\(t)  (2-3) 

With  the  equations  of  motion  expressed  in  this  form,  a  state  space  rep¬ 
resentation  can  easily  be  attained  (Ref  15)  and  therefore  the  system  can 
be  modeled  as 


X(t)  -  lA]X(t)  +  (B]u(t) 
T(t)  -  (C]X(t) 


(2-4) 

(2-5) 


. 


■a  ■  ^ 


where  the  state  matrices,  state  vector  and  control  vector  are  defined  as 
in  Appendix  A. 


Optimal  Regulator  Theory 

With  the  aeroelastic  equations  of  motion  modeled  as  in  Eq  2-4,  the 
control  law  synthesis  procedure  can  be  treated  as  a  linear  optimal  reg¬ 
ulator  problem.  The  function  to  be  minimized  is  an  infinite  time  integral, 
quadratic  cost  function  in  terms  of  the  states  and  controls: 


a 


j 


[XT(t)[Q]X(t)  +  uT(t)(R]u(t)]dt 


(2-6) 


where  [Q]  and  [R]  are  weighting  matrices  on  the  states  and  controls  re¬ 
spectively.  [Q]  is  chosen  to  be  positive  semi-deflnlte  and  [R]  chosen 
to  be  positive  definite,  which  guarantees  a  stable  feedback  control  law. 
The  minimization  of  this  performance  index  leads  to  the  optimal  control 
law 


u*(t)  -  lG]X(t) 


(2-7) 


where 


[G]  -  -IR]-1  IBJ T  [P] 


(2-8) 
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For  the  time-invariant  case  (constant  coefficient  differential  equations), 
[P]  is  the  steady  state  solution  to  the  Matrix  Ricatti  equation: 


-IP]  -  IP]  [a]  +  U1T[P1  -  [PHB][R]*1(B]T(PJ  +  [Q]  -  (OJ  (2-9) 

The  application  of  quadratic  optimization  is  an  Iterative  procedure  of 
selecting  the  appropriate  performance  function  through  changes  in  the 
weighting  matrices  [Q]  and  [&] .  The  choice  of  the  approplate  weighting 
matrices  is  dependent  on  the  designer's  past  experience  and  his  under¬ 
standing  of  the  physics  of  the  problem. 


Control  Law  Synthesis 

Controller  Design.  To  begin  the  optimal  control  law  synthesis  pro¬ 
cedure,  it  was  desired  to  uncouple  the  states  of  the  system  by  utilizing 
a  state  transformation  to  diagonalize  the  state  coefficient  matrix  [A]. 
Employing  the  substitution, 

X(t)  -  [T]Z(t)  (2-10) 

where  the  columns  of  {T]  are  the  eigenvectors  of  [A]  (Ref  16),  the  orig¬ 
inal  system  represented  by  Eqs  2-4  and  2-5  can  be  transformed  to: 

(T]Z(t)  -  tA](TjZ(t)  +  (Bju(t) 

Y(t)  -  (CHT)Z(t) 


or 
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I  A]2(t)  +  lB']u(t, 
lC'JZ(t) 


(2-11) 


Z(t)  “ 
I(t)  - 


where 

l  A]  -  iTf1  [A] It] 

ib» j  -  nr1  is] 

end 

1C)  -  IC][T] 


a 


With  the  system  modeled  In  this  form,  en  author-modified  version  of  the 
program  OPT CON  (Ref  20)  was  employed  to  solve  the  matrix  Rleattl  equation 
given  by  Eq  2-9,  and  obtain  the  optimal  control  law  u*(t)  ■  [G]X(t)  . 

Due  to  the  large  range  of  numerical  values  spanned  when  the  controllabil¬ 
ity  matrix  is  formed  Inside  0PTC0N,  the  Ricattl  solver  subroutine  could 
not  function  properly.  This  resulted  in  the  program  defining  the  system 
as  being  uncontrollable  when  this  was  clearly  not  the  case. 

To  eliminate  the  Inherent  numerical  difficulties,  it  was  decided 
that  the  three  very  high  frequency  aerodynamic  lag  states  would  be  neg¬ 
lected  for  the  feedback  gain  matrix  calculation.  In  general,  this  would 
not  be  necessary  if  a  more  efficient  Ricattl  equation  solver  could  be 
utilised.  The  system  represented  by  Sq  2-11  can  now  be  presented  as: 


2  (t) 

A.  0  " 

2  (t) 

B,  ' 

—  1 

m 

1 

“1 

+  * 

1 

2  (t) 

0  A 

2  (t) 

* 

B  ' 

-2 

L  2j 

-2 

2 

u(t) 


(2-12) 


where  A 2  is  a  3x3  matrix  containing  the  high  frequency  roots.  The  control 
law  was  then  synthesized  based  on  the  following  reduced  order  subsystem 
of  Eq  2-12: 


at 


i2(t)  -  U2]Z2(t)  +  lB2’]u(t)  (2-13) 

The  above  IA2)  and  I B matrices  were  then  input  into  the  modified  OPTCON 
program  and  the  optimal  control  law 

u'(t)*  -  lG2]Z2<t)  (2-14) 

was  determined  where  [G  ]  Is  formed  as  in  Eq  2-8.  This  control  law  was 

2 

then  Incorporated  into  Eq  2-11  by  letting: 

u*(t)  -  [6  )Z(t)  (2-15) 


where 


[G  ]  •  [  0  :  G  J 
1  2 


(2-16) 


Mow  the  closed  loop  system  takes  on  the  form: 


Z(t)  -  [A  )Z(t)  +  [B']u*(t) 

-  t A ]Z(t)  +  IB* ) [G1 ]Z(t) 

-  <[A]  +  IB'JIG  ])Z(t) 

-  [A’  JZ(t)  (2-17) 

CL  ~ 


11 


Ol 


The  closed  loop  eigenvalues  are  then  determined  by  the  matrix  [A'  ]. 

vL 

The  optimal  control  law  given  by  Eq  2-15  can  also  be  put  back  into 
the  original  system  by  using  the  inverse  relationship; 


Z(t)  -  [T]_1X(t) 


(2-18) 


resulting  in: 


u*(t)  -  IGJITI^XCO 
-  (GjX(t) 


(2-19) 


Therefore, 

X(t)  -  [A]X(t)  +  [Bju* (t) 

-  (A]X(t)  +  [B] [G]X(t) 

-  ([A]  +  [B][G])X(t) 

-  [A^JXU)  (2-20) 

and  the  resulting  eigenvalues  of  l ]  and  [A^]  are  nearly  identical. 

The  coupling  of  the  states  resulting  from  the  Inverse  [Tj  transformation 
along  with  the  numerical  accuracy  of  the  operations  involved  does  cause 
some  original  system  closed  loop  eigenvalues  to  shift  slightly. 

The  optimal  control  law  given  by  Eq  2-15  or  Eq  2-19  assumes  that  it 
is  possible  to  feed  back  the  complete  state  vector.  In  actuality,  this 
Is  not  possible  as  the  number  of  states  that  are  actually  measurable  are 
limited  by  the  number  of  sensors  used  and  how  sophisticated  these  sensors 
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are.  For  this  study,  two  Independent  sensors  were  utilized;  one  to  sense 
wing  vertical  displacement  and  one  to  sense  wing  tip  angular  acceleration. 
The  results  of  Moll's  analyses  indicated  that  a  leading  edge  surface  com¬ 
manded  by  displacement  feedback  provided  a  reasonable  control  system  de¬ 
sign  for  preventing  divergence  of  the  cantilever  wing  or  the  body  freedom 
flutter  instability  associated  with  the  model  free  in  pitch.  The  dis¬ 
placement  sensor  was  positioned  near  the  intersection  of  the  wing  second 
bending  node  line  and  the  wing  torsion  node  line.  This  location  was  de¬ 
termined  to  be  optimum  for  feeding  back  the  bending  motion  of  the  first 
elastic  mode  with  minimum  inputs  from  all  other  elastic  modes.  Analyses 
also  indicated  that  a  trailing  edge  control  surface  commanded  by  angular 
acceleration  of  the  wing  tip  perpendicular  to  the  elastic  axis  provided 
an  adequate  input  for  controlling  the  bending/ torsion  flutter  mode. 

Feeding  back  wing  tip  acceleration  assured  maximum  input  from  the  torsion 
mode  with  minimum  response  from  the  bending  modes. 

Since  the  full  state  vector  is  not  available  through  direct  measure¬ 
ment,  it  is  necessary  to  design  a  deterministic  observer  (state  estimator) 
to  reconstruct  the  state  vector  from  the  limited  amount  of  sensor  infor¬ 
mation  available. 

Observer  Design.  In  order  to  reconstruct  the  state  X(t)  of  the  or¬ 
iginal  system  from  the  observed  output  vector  Y(t)  as  given  by  Eq  2-5,  a 
linear  differential  system  must  be  formed  whose  output  is  to  be  an  approx¬ 
imation  to  the  state  X(t).  It  will  be  Illustrated  in  the  following  devel- 
opement  as  to  what  structure  this  system  should  have  and  how  it  should 


behave 


at 


The  n-dimensional  system  represented  by 

X(t)  -  [F)X(t)  +  [G]Y(t)  +  [H]u(t) 


where 


(2-21) 


[FJ  -  [A]  -  IK]  [C] 

[GJ  -  IK] 

IH]  -  IB]  (2-22) 

is  a  full-order  observer  for  the  n-dimensional  system  given  in  Eq  2-4  if 

lv  -  s<*0> 

implies 

X(t)  -  X(t)  t  i  tQ  , 

for  all  u(t),  t  2  tg  (Ref  11). 

The  observer  given  by  Eq  2-21  is  called  a  full-order  observer  since 
its  state  X(t)  has  the  same  dimension  as  the  state  X(t)  of  the  original 
system.  With  the  system  matrices  defined  as  in  Eq  2-22,  the  observer 
structure  becomes: 

|(t)  -  UJXU)  +  tB]u(t)  +  lK](Y(t)  -  (C]X(t)) 
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for  all  e^t^),  If  and  only  if,  Che  observer  is  asymptotically  stable. 
Comparing  Eq  2-23  and  Eq  2-25,  we  see  that  the  stability  of  the  observer 
and  the  asymptotic  behavior  of  the  reconstruction  error  are  both  deter¬ 
mined  by  the  behavior  of  the  matrix  ((A]  -  [K][C]).  This  clearly  shows 
that  the  reconstruction  error  ,e(t)  approaches  zero,  independent  of  its 
initial  value,  if  and  only  if  the  observer  is  asymptotically  stable. 

This  is  a  very  desirable  characteristic  of  a  properly  designed  observer. 

The  preceding  developement  clearly  illustrates  that  deterministic 
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observer  design  revolves  about  determining  the  gain  matrix  [K]  such  that 
the  reconstruction  error  differential  equation  is  asymptotically  stable. 
Since  in  this  analysis  we  are  dealing  with  time-invariant  coefficient 
matrices,  the  stability  of  the  observer  is  determined  by  the  location  of 
the  eigenvalues  of  the  matrix  ((A]  -  [KltC]).  We  refer  to  these  eigen¬ 
values  as  the  observer  poles.  It  can  also  be  proven  (Ref  11)  that  all 
observer  poles  can  be  arbitrarily  located  in  the  complex  plane  by  a  suit¬ 
able  choice  for  (K]  if  the  system  is  completely  recons tractable.  In  gen¬ 
eral,  it  is  desirable  to  choose  [X]  such  that  the  observer  poles  are  deep 
in  the  left-half  complex  plane  so  as  to  obtain  fast  convergence  of  the 
reconstruction  error  to  zero.  This  generally  requires  the  gain  matrix 
[K]  to  be  large  which  could  cause  the  observer  to  be  very  sensitive  to 
observation  noise  that  may  be  present. 

With  the  reconstruction  error  defined  as  in  Eq  2-24,  the  optimal  con¬ 
trol  law  based  on  the  estimate  of  X(t)  can  be  represented  as: 

u*(t)  -  [GjX(t) 

-  [G](X(t)  -  e(t))  (2-26) 

When  this  control  law  is  Incorporated  back  into  the  original  system,  the 
following  state  equation  is  obtained: 

X(t)  -  !AlX(t)  +  (Bju* (t) 

-  (AjX(t)  +  lB]{G](X(t)  -  e(t)) 

-  UA]  +  [B)[G})X(t)  -  [B}[G]e(t)  ; 

-  |ACL)X(t)  -  lB](G]e(t) 


(2-27) 


As  previously  discussed,  for  a  properly  designed  observer  the  error  £(t) 
approaches  zero  over  time  and  the  original  state  system  would  be  attained. 

As  with  the  controller,  the  observer  design  can  be  based  on  the  un- 
coupled  Z(t)  system  for  computational  ease.  A  vector  £(t)  is  defined  to 
be  an  estimate  of  £(t)  such  that 

Z(t)  -  lAl£(t)  +  [B*]u(t)  +  (K  )(Y(t)  -  Y(t)) 

Y(t)  -  IC’JZU)  (2-28) 

A  new  reconstruction  error  for  the  diagonalized  system  is  defined  to  be: 

e'(t)  -  Z(t)  -  Z(t)  (2-29) 

Then  by  taking  the  difference  between  Eq  2-11  and  Eq  2-28  the  resulting 
state  model  for  the  reconstruction  error  can  be  represented  as 

e’(t)  -  [  A]e'(t)  -  iKjHC'Je’U) 

-  (I  A  ]  -  l^JtC'DeUt)  (2-30) 

and  the  stability  of  the  error  is  determined  by  the  eigenvalues  of  the 
([A]  -  lK^}[C'))  matrix.  Since  the  objective  is  to  determine  [K^]  such 
that  the  reconstruction  error  is  asymptotically  stable,  existing  software 
can  be  utilized  by  forming  an  auxiliary  problem  using  the  duality  relation¬ 
ship  of  Eq  2-28: 

i’(t)  -  lA]Ti’(t)  +  [C'J  M(t)  (2-31) 
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The  control  can  be  determined  using  the  techniques  previously  discussed 
to  obtain  the  desired  eigenvalue*  of  the  reconstruction  error,  and  hence 
the  observer.  Therefore,  M(t)  can  be  represented  as 


£#(t)  -  £'(t) 


(2-32) 


and  the  eigenvalues  of  the  closed  loop  reconstruction  error  system  are 

T  T 

determined  by  the  ([ A ]  -  (C'J  [K^J  )  matrix. 

Once  again,  to  prevent  numerical  difficulties  inside  the  modified 
OPT CON  program,  it  was  oecessary  to  obtain  the  observer  control  law  using 
a  reduced  order  model  and  the  same  methodology  utilized  in  the  controller 
design.  Therefore,  the  following  reduced  order  auxiliary  problem  is 
formed: 


JL'(t)  -  t  A  ]£’2<t)  +  tc']  V(t) 


and  the  control  law  is  given  by 


(2-33) 


it  (t)  -  -[k2j  i2(t) 


(2-34) 


The  gain  matrix  [R^]  utilized  in  Eqs  2-30  and  2-32  to  determine  the  ob¬ 
server  poles  can  be  formed  from  the  l^)  gain  matrix  in  a  similar  manner 
that  the  controller  control  law  was  formed  in  Eq  2-16: 


T  T 

1^1  J 
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ux]  - 


(2-35) 


Once  this  gain  matrix  is  determined,  the  reconstruction  error  can  be  ex¬ 
pressed  in  terms  of  the  original  non-diagonalized  system  by  using  the  re¬ 
lationship  given  in  Eq  2-10,  that  is 


Z(t)  -  lT]_1X(t)  and  Z(t)  -  [T]_1X(t) 


(2-36) 


e'(t)  -  Z(t)  -  i{t) 

“  [T] _1(X(t)  -  X(t)) 
-  lTl”1e(t) 


(2-37) 


Therefore,  Eq  2-30  can  be  rewritten  as: 


i'(t)  - 
It]  _1e(t)  - 


(U  ]  -  [KjHC’De'U) 


le(t)  -  (1  A  ]  -  UjllC'DlTj-^t) 
e(t)  -  IT] (| A ]  -  [K1]tC'])(T]"1e(t) 
*  Ue]e(t) 


(2-38) 


Finally,  an  augmented  state  vector  for  the  entire  plant-observer  system 
can  be  formed  using  Eqs  2-27  and  2-38: 


X(t)  f  ArT  -BG|  X(t) 


Z  (t)  - 
CL  e(t) 


(2-39) 


Ae  e(t) 
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This  augmented  state  model  clearly  illustrates  that  the  closed  loop  sys- 
tern  poles  consist  of  two  distinct  sets;  the  poles  of  ([A]  +  [B)[G])  which 
determine  the  behavior  of  the  state  X(t),  and  the  poles  of  ([A]  -  (K][C]> 
which  determine  the  behavior  of  the  error  e_(t).  This  is  known  as  the  Sep- 
aration  Theorem  (Ref  14:544)  and  illustrates  that  the  determination  of  a 
feedback  control  law  and  the  construction  of  an  observer  can  be  considered 
as  separate  problems. 


III.  Results 


As  outlined  In  Appendix  A,  the  elements  of  the  state  coefficient 
matrices  are  functions  of  velocity  and  density,  and  therefore  for  a  fixed 
altitude  the  characteristic  roots  can  be  determined  as  velocity  is  in¬ 
creased  from  zero.  The  velocity  at  which  the  real  part  of  one  of  the 
roots  becomes  zero  is  the  divergence  or  flutter  velocity  depending  on 
whether  the  frequency  Is  zero  or  nonzero  when  it  crosses  the  Imaginary 
axis.  Presented  in  Figure  2  is  a  velocity  root  locus  for  the  open  loop, 
unaugmented  system  at  sea  level.  It  was  decided  to  present  this  plot  as 
a  function  of  equivalent  velocity  so  that  a  similar  plot  can  be  made  for 
any  desired  altitude  by  multiplying  the  velocity  by  the  square  root  of 
the  density  ratio.  It  should  be  mentioned  here  that  although  these  plots 
are  called  "root  locus  plots,"  they  do  not  adhere  to  traditional  root 
locus  construction  rules  and  Just  serve  to  illustrate  the  movement  of  the 
characteristic  roots  as  velocity  is  varied.  The  passive  divergence  speed 
(Vq)  was  predicted  to  occur  at  115  fps  and  the  passive  bendiz^/ torsion 
flutter  speed  (V^)  was  predicted  to  occur  at  156  fps  at  a  flutter  freq¬ 
uency  (w^)  of  104.9  rad/sec.  For  the  augmented  analyses,  the  goal  was  to 
determine  an  optimal  control  law  for  the  design  flight  condition  which 
increased  the  onset  of  instability  20Z  above  the  unaugmented  wing  bend¬ 
ing/torsion  flutter  speed.  To  accomplish  this  goal,  the  state  coefficient 
matrices  were  calculated  at  the  design  velocity  of  1.2V^  “  167  fps  and 
utilized  in  the  control  law  synthesis.  Wit)-  the  state  model  formed,  the 
open  loop  eigenvalues  at  this  design  condition  were  calculated  and  are 
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contained  In  Table  1 


Table  1 

Open  Loop  Eigenvalues  at  Design  Velocity 


1  21.190  0.0 

2  -57.342  0.0 

3  -14.034  74.807 

4  -14.034  -74.807 

5  10.393  87.113 

6  10.393  -87.113 

7  -154.03  0.0 

8  -178.77  0.0 

9  -179.75  0.0 

10  -76714.  0.0 

11  -76969.  0.0 

12  -78118.  0.0 


Ol 

Throughout  the  control  law  synthesis  procedure,  the  control  weight¬ 
ing  matrix  [R]  was  set  equal  to  the  identity  matrix  and  only  the  state 
weighting  matrix  lQ]  was  varied  to  obtain  the  desired  closed  loop  behav¬ 
ior.  This  is  an  acceptable  approach  as  it  is  the  relative  relationship 
between  [Q]  and  [R]  that  Is  of  Importance.  Initially,  the  state  weight¬ 
ing  matrix  [Q]  was  set  equal  to  zero  since  this  yields  a  set  of  gains 
that  are  "cheapest"  in  terms  of  input  amplitude  (Ref  11:  289).  Table  2 
contains  the  closed  loop  eigenvalues  at  the  design  condition  for  [Q] 
equal  to  zero.  It  can  be  seen  that  the  effect  of  this  control  law  is 
to  leave  all  stable  eigenvalues  unchanged  and  reflect  the  unstable  eigen¬ 
values  about  the  imaginary  axis.  It  should.' also  be  noted  that  this  con¬ 
trol  law  assumes  that  the  full  state  vector  is  available  for  feedback. 
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Table  2 


Closed  Loop  Eigenvalues  Assuming  Full  State 
Feedback  Gains  at  Design  Velocity 


Mo 

Real 

Imaginary 

1 

-21.190 

0.0 

2 

-57.342 

0.0 

3 

-14.034 

74.807 

4 

-14.034 

-74.807 

5 

-10.393 

87.113 

6 

-10.393 

-87.113 

7 

-154.03 

0.0 

8 

-178.77 

0.0 

9 

-179.75 

0.0 

10 

-76714. 

0.0 

11 

-76969. 

0.0 

12 

-78118. 

0.0 

Table  3 

Closed  Loop  Eigenvalues  at  Off-design  Points 


Eigenvalues 

Velocity(fps) 

Real 

Imaginary 

94  <.6Vf) 

-2.843 

12.834 

-2.843 

-12.834 

-2.283 

78.904 

-2.283 

-78.904 

-1.290 

135.550 

-1.290 

-135.550 

-85.303 

0.0 

-90.023 

0.0 

-90.381 

0.0 

-38562. 

0.0 

-38690. 

0.0 

-39268. 

0.0 

115  (VD) 

-279.79 

0.0 

-9.607’ 

0.0 
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Table  3  (Coat'd.) 


Velocity(fps) 

Real 

Eigenvalues 

Imaginary 

115  (coat'd.) 

-3.026 

79.076 

-3.026 

-79.076 

-1.707 

127.740 

-1.707 

-127.740 

-102.58 

0.0 

-110.06 

0.0 

-110.56 

0.0 

-47177. 

0.0 

-47333. 

0.0 

-48040. 

0.0 

156  (V,) 

-31.002 

0.0 

I 

-30.813 

0.0 

-6.296 

80.479 

-6.296 

-80.479 

-0.200 

108.330 

-0.200 

-108.330 

-133.63 

0.0 

-149.19 

0.0 

-149.96 

0.0 

-63997. 

0.0 

-64212. 

0.0 

-65167. 

0.0 

203  (1.3V,) 

-79.674 

0.0 

£ 

-162.99 

0.0 

10.503 

72.841 

10.503 

-72.841 

-13.000 

68.953 

-13.000 

-68.953 

-194.05 

0.0 

-195.13 

0.0 

-392.86 

0.0 

-83280. 

0.0 

-83556. 

0.0 

-84819. 

0.0 
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Table  3  (Cont'd.) 


I  _ 

Eigenvalues 


Velocity(fps) _ Real _ Imaginary 


219  (1.4V  )  -106.17  0.0 

f  -106.58  0.0 

I  -9.270  65.287 

-9.270  -65.287 

-6.175  61.454 

-6.175  -61.454 

-170.08  0.0 

-209.3:  0.0 

|  -210.51  0.0 

-89841.  0.0 

-90155.  0.0 

-91513.  0.0 


I 


I  a 

i 
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To  assess  the  robustness  of  the  feedback  control  law,  the  feedback 
gain  matrix  calculated  at  the  design  flight  condition  with  [Q]  equal  to 
the  zero  matrix  was  applied  at  off-design  flight  conditions  assuming 
full  state  feedback.  Contained  in  Table  3  are  the  dosed  loop  eigen¬ 
values  at  representative  off-design  velocities.  As  shown  by  this  table 
and  in  Figure  3,  this  gain  matrix  does  perform  adequately  up  to  the  de¬ 
sign  velocity  of  187  fps  and  slightly  beyond,  but  narrowly  maintains 
stability  at  156  fps.  Above  approximately  194  fps,  the  torsion  mode 
becomes  unstable  but  does  recross  the  imaginary  axis  into  the  stable 
regime  at  higher  velocities.  This  is  not  an  unusual  result  as  the  con¬ 
trol  law  is  optimized  for  one  flight  condition  and  does  not  guarantee 
acceptable  performance  at  any  other  flight  condition.  The  fact  that  the 
closed  loop  system  is  stable  up  to  the  design  velocity  is  an  excellent 
result,  but  to  achieve  robustness  throughout  the  complete  flight  regime 
this  control  would  not  be  acceptable. 


1 

.( 

1 
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An  attempt  was  then  made  to  vary  [ Q ]  and  the  individual  state  weight¬ 
ings  inside  [Q]  to  obtain  &  feedback  gain  matrix  that  exhibits  satisfact¬ 
ory  performance  for  velocities  greater  than  187  fps.  Since  it  is  the 
torsion  mode  that  becomes  unstable  at  higher  velocities ,  the  state  assoc¬ 
iated  with  that  mode  was  weighted  more  heavily  inside  the  weighting  mat¬ 
rix  [Q] .  This  has  the  effect  of  directing  more  of  the  available  control 
authority  to  controlling  this  mode  of  instability.  After  several  iter¬ 
ations,  a  weighting  of  100  on  the  torsion  mode  and  0.1  on  the  other  modes 
produced  more  desirable  results  as  shown  in  Figure  4.  This  feedback  con* 
trol  law  exhibited  satisfactory  performance  for  all  of  the  velocities  in¬ 
vestigated. 

It  was  not  the  intent  of  this  thesis  to  assess  whether  this  closed 
loop  performance  is  acceptable  when  considering  the  actual  dynamic  re¬ 
sponse  of  the  aircraft  when  disturbed  from  a  steady  state  condition.  It 
should  be  understood  that  this  study  involved  standard  eigenvalue  anal¬ 
ysis  techniques,  and  the  criterion  to  judge  acceptable  performance  was 
only  that  the  eigenvalues  have  negative  real  parts,  i.e.  the  system, 
when  excited  from  a  steady  state  condition,  would  have  a  transient  re¬ 
sponse  which  decays  to  zero  over  time. 

As  discussed  earlier,  the  full  state  vector  is  not  available  for 
measurement  and  therefore  an  observer  (state  estimator)  must  be  designed 
to  provide  an  estimate  of  the  state  vector  to  be  used  in  the  feedback 
control  law.  The  observer  gains  were  determined  using  the  diagonalized 
state  coefficient  matrix  and  the  transpose  of  the  output  coefficient 
matrix  which  were  input  into  the  modified  0PTC0N  program,  and  the  weight¬ 
ing  matrices  varied  to  obtain  the  desired  observer  poles.  Once  again, 
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Figure  4.  Closed  Loop  Velocity  Root  Locus  With  Improved  Weightings 
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Che  control  weighting  matrix  { 0B  was  kept  equal  to  the  identity  matrix 
and  only  [Q]qb  was  varied.  Alter  several  iterations,  the  largest  weight¬ 
ing  matrix  that  could  be  utilized  was  [Q]  “  [I]  due  to  numerical  dif- 

OB 

ficulties  Inside  the  OPTCON  program.  It  was  decided  that  although  the 
observer  poles  were  not  as  deep  in  the  left  half  plane  as  was  desired, 
the  eigenvalues  obtained  were  adequate  for  the  current  analysis  and  for 
demonstration  of  optimal  control  theory  methooology. 

The  observer  poles  obtained  for  this  case  at  the  design  velocity 
are  presented  in  Table  4  and  can  be  seer  to  be  adequately  stable,  but  not 


Table  4 

Observer  Poles  for  Design  Velocity 
"  ^OB  “  ^ 


No 

Seal 

Imaginary 

1 

-25.257 

0.0 

2 

-62.662 

0.0 

3 

-14.943 

75.729 

4 

-14.943 

-75.729 

5 

-10.943 

87.413 

6 

-10.943 

-87.413 

7 

-158.52 

0.0 

8 

-188.38 

0.0 

9 

-3846.3 

0.0 

10 

-76714. 

0.0 

11 

-76969. 

0.0 

12 

-78118. 

0.0 

as  rapidly  convergent  as  was  desired.  As  with  the  closed  loop  eigenvalues, 
these  observer  poles  would  vary  as  velocity  changes,  and  improved  weight¬ 
ing  matrices  might  be  required.  This  was  not  Investigated  in  this  study 
due  to  the  numerical  problems  encountered  in  the  observer  design.  In  the 
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procedure  to  design  an  observer  to  be  actually  implemented,  it  might  be 
warranted  to  schedule  the  observer  gain  matrix  as  a  function  of  dynamic 
pressure  so  as  to  assure  fast  convergence  of  the  reconstruction  error  to 
zero  throughout  the  flight  envelope. 

One  other  area  of  Interest  was  directed  at  determining  the  effect  of 
neglecting  the  three  high  frequency  aerodynamic  states  in  the  control  law 
synthesis.  This  results  in  a  reduced  order  state  model  (9th  instead  of 
12th)  which  simplifies  the  numerical  calculations  and  improves  the  math¬ 
ematical  Inaccuracies.  When  the  state  model  was  formed,  the  feedback 
gain  matrix  calculated  and  the  closed  loop  eigenvalues  determined,  it 
was  discovered  that  neglecting  these  poles  had  no  noticable  effect  on 
the  closed  loop  performance  (see  Table  5).  This  is  as  would  be  expected 


Table  5 

Reduced  Order  Model,  Closed  Loop  Eigenvalues 
At  Design  Velocity 


No 

Real 

Imaginary 

1 

-21.190 

0.0 

2 

-57.342 

0.0 

3 

-14.034 

74.807 

4 

-14.034 

-74.807 

5 

-10.393 

87.113 

6 

-10.393 

-87.113 

7 

-154.03 

0.0 

8 

-178.77 

0.0 

9 

-179.75 

0.0 

due  to  the  large  magnitude  of  the  neglected  poles. 

This  is  also  a  very  Interesting  result  as  it  illustrates  that  so- 


phisticated,  three  dimension^.-  unsteady  aerodynamic  modeling  techniques 
can  be  utilized  to  form  the  aeroelastic  equations  of  motion,  but  sub¬ 
sequently  the  equations  associated  with  the  high  frequency  aerodynamic 
lag  states  can  be  neglected.  This  has  the  effect  of  reducing  the  order 
of  the  state  model  formed  and  thereby  reducing  the  complexity  of  the 
numerical  operations  performed,  and  eliminates  possible  numerical  dif¬ 
ficulties  inside  the  Rlcatti  equation  solver  subroutine. 


IV.  Conclusions  and  Recommendations 


It  has  been  demonstrated  herein  that  optimal  control  theory  tech- 
nlques  can  be  applied  adequately  for  the  prevention  of  aeroelastic  insta¬ 
bilities  on  a  forward  swept  wing  aircraft.  The  feedback  control  law  syn¬ 
thesised  at  the  chosen  design  point  performed  adequately  throughout  the 
flight  regime.  It  was  also  illustrated  that  an  observer  could  be  proper¬ 
ly  designed  to  reconstruct  the  original  state  vector  from  available  sensor 
measurements.  The  feedback  control  law  was  also  formed  from  a  reduced 
order  state  model  which  neglected  the  equations  associated  with  the  high 
frequency  aerodynamic  lag  states.  The  control  law  synthesized  from  this 
state  model  performed  similarly  for  all  velocities  investigated,  and  the 
f S  resulting  closed  loop  eigenvalues  were  not  altered  by  any  noticeable 

amount. 

The  results  of  this  study  are  just  a  foundation  for  an  area  of  con¬ 
trol  system  design  that  warrants  further  investigation.  It  is  planned  in 
in  the  near  future  to  expand  this  study  to  Include  similar  analyses  for 
the  forward  swept  wing  model  free  in  pitch  and  subsequently  for  the  model 
free  in  pitch  and  plunge.  For  the  model  free  in  pitch,  coupling  between 
the  wing  first  bending  mode  and  the  rigid  body  pitch  mode  occurs  at  a  vel¬ 
ocity  lower  than  the  static  divergence  speed,  and  therefore  becomes  the 
instability  of  interest.  To  complete  the  analysis  of  the  forward  swept 
wing  model,  the  case  for  the  model  free  in  pitch  and  plunge  will  be  in¬ 
vestigated  so  as  to  actively  control  the  more  realistic  representation 
of  an  actual  aircraft  employing  a  forward  swept  wing  design. 
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There  ere  several  additional  concepts  that  could  be  applied  to  this 
study  to  further  enhance  the  validity  of  the  results.  Although  it  is 
realized  that  control  surface  displacements  and  rates  are  important  in 
the  design  and  effectiveness  of  a  control  law,  the  calculation  of  these 
quantities  was  not  considered  in  this  study.  Also,  sensor  and  actuator 
dynamics  were  neglected  in  this  analysis.  The  control  of  divergence  for 
the  cantilever  wing  and  the  body  freedom  flutter  instability  for  the  model 
free  in  pitch  would  be  expected  to  be  Insensitive  to  actuator  and  sensor 
dynamics  because  of  the  low  frequencies  involved.  However,  It  is  antic** 
ipated  that  the  control  of  the  higher  frequency  wing  bending/torsion  flut¬ 
ter  mode  would  be  affected  by  the  addition  of  these  components. 

Another  possible  area  of  investigation  slight  be  to  perform  a  gust 
analysis  on  the  closed  loop  control  system  to  obtain  a  more  robust  con¬ 
trol  law  which  is  adequate  throughout  the  flight  envelope.  Several  tech¬ 
niques  have  been  investigated  such  as  injecting  white  noise  into  the  con¬ 
troller  system  model  at  the  point  of  entry  of  the  control  inputs,  during 
the  process  of  tuning  the  filter  (observer).  This  would  allow  the  con¬ 
troller  to  exhibit  more  robustness  at  the  design  points,  which  may  reduce 
the  number  of  required  design  conditions  necessary  to  have  adequate  per¬ 
formance  throughout  the  flight  envelope.  It  may  even  be  desirable  to 
design  a  fixed  controller  which  Is  sufficiently  robust  to  produce  accept¬ 
able  performance  for  all  conditions  in  the  aircraft's  flight  envelope. 

A  natural  extension  to  this  is  to  consider  injecting  time-correlated  noise. 
This  allows  the  robustlflcation  technique  to  be  applied  only  over  a  de¬ 
sired  frequency  range  rather  than  over  all  frequencies  as  in  the  previous 
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Appendix  A:  Aeroelaatlc  Equations  of  Motion 


The  aeroelastic  equations  of  notion  of  a  flexible  aircraft  in  an 
airstream  can  be  represented  as  (Ref  12): 

lM]^(t)  +  [C]^(t)  +  [K]£(t)  -  F(t)  (A-l) 

where  [M] ,  (C]  and  [K]  are  the  generalized  mass,  damping  and  stiff nesa 
matrices  obtained  using  a  set  of  generalized  coordinates  ^(t)  and  sev¬ 
eral  natural  vibration  mode  shapes,  and  £(t)  represents  the  unsteady 
aerodynamic  forces.  In  Noll's  study  (Ref  9),  the  forces  were  obtained 
from  subsonic  doublet  lattice  unsteady  aerodynamic  theory  (Ref  10)  and 
are  defined  to  be 

F(t)  -  -l/2pV2s[Q(k)]£(t)  (A-2) 

where  a  represents  a  reference  area. 

The  generalized  aerodynamic  force  coefficient  matrix  [Q(k)]  was 
computed  from 


AP.  h- 

- %  — 

l/2pV2  s 


(A-3) 


where  h^  is  the  normalized  vertical  displacement  in  the  1th  mode.  The 
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coefficient  represents  the  non-dimensional  modal  force  in  the  ith  mode 
due  to  pressure  from  the  jth  mode  and  is  dependent  on  the  reduced  freq¬ 
uency  k  (where  k  *  bw/V).  The  expression  for  F(t)  can  now  be  substituted 
Into  Eq  A-l  and  the  Laplace  Transform  taken  to  yield  (for  zero  initial 
conditions):  • 


( [M]s2  +  [C]s  +  IK]  +  l/2pV2s[Q(s)])£(s)  -  0 


(A-4) 


Before  the  equations  of  motion  can  be  transformed  into  the  time  domain. 
It  is  necessary  to  express  the  elements  of  the  [Q(s>]  matrix  in  a  con¬ 
venient  form.  For  the  current  analysis,  it  was  assumed  that  the  ele¬ 
ments  of  [Q(8>]  are  given  by 


S  -  bs/V 

which  is  a  Pade'  approximant  representation  of  the  aerodynamic  forces 
(Ref  22).  The  coefficients  and  were  obtained  using  a  least  squares 
fleeing  process  over  a  reduced  frequency  range  of  ineerese.  The  gen¬ 
eralized  force  matrix  now  becomes  (in  terms  of  the  Laplace  variable  s): 

(Cft]V3  +  lC,]V2bs  +  lC,]Vb2s2  +  1C ,lb3s3 

|Q(8>]  -  - 2 - r-J - = - 2 - - - 1 -  (A-6) 

VJ  +  DjbV^s  +  DjVb^s 
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i  a 


Substituting  this  expression  into  Eq  A-4  and  multiplying  through  by  the 
denominator  polynomial,  the  transformed  equations  of  nfbtion  take  the 
form: 


<U4]s4  +  U3]s3  +  [A2]s2  +  [A^s  +  lA0])i(s)  *  0 


(A-7) 


The  elements  of  the  (A^J  matrices  are  functions  of  velocity  and  there¬ 
fore  the  roots  of  the  characteristic  equation  can  now  be  determined 
as  velocity  is  varied. 

When  control  surfaces  are  added  to  the  flexible  wing,  their  effect 
must  be  taken  into  consideration  when  forming  the  equations  of  motion. 

We  can  begin  the  derivation  by  incorporating  the  effects  of  control  sur¬ 
faces  into  Eq  A-l  : 


[M)£(t)  +  [Cj£(t)  +  lK]A(t)  +  l/2pV2s[Q(k)]£(t) 


+  !Mc]j^,(t)  +  l/2PV2slQc(k)]£c(t)  -  0 


(A-8) 


where  [M  ]  and  [Q  ]  are  matrices  of  order  NM  x  NC  (NM  equals  the  number 
c  c 

of  elastic  modes  and  NC  is  the  number  of  control  surfaces).  For  this 
study,  two  active  control  surfaces  were  employed  (leading  and  trailing 
edge)  resulting  in 


S.CM  -  UTE<t)  sLE(t)I 


(A-9) 


The  (Qc(s)]  matrix  to  be  used  in  the  Laplace  Transform  of  Eq  A-8  is  also 
obtained  using  a  Pade’  approximant  for  the  control  surface  aerodynamics 
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and  takes  on  the  same  form  as  Eq  A- 6  : 


lCn  ]V3  +  [C.  ]V2bs  +  [c  jvb2s2  +  [C,  JbV 
IQ  (s)1  -  — - - 2fi - Is - 

l  Xc  \  ®/  J  o  ry  2  2 

V3  +  DjbVzs  +  D2Vb  s 


(A-10) 


The  denominator  of  the  lQc(s)]  matrix  is  also  constrained  to  be  the  same 
as  the  denominator  of  the  [Q(s)J  matrix. 

Mow,  when  the  aerodynamic  modeling  is  substituted  into  the  Transform 
of  Eq  A-8  and  the  resulting  equation  is  multiplied  through  by  the  denom¬ 
inator  of  Eqs  A-6  and  A-10,  the  aeroelastic  equations  of  motion  take  on 
the  form: 


<[AaJs4  +  lA^s3  +  [A^s2  +  [A^a  +  [AQJ)£(s) 

+  (IBaJs4  +  [B31s3  +  [B2Js2  +  [Bjs  +  tB()l>icCs)  -  0  (A-ll) 


where 

|Aa)  -  [MlD2b2 

U3J  -  [M]DjbV  +  [CjD2b2  +  l/2pslC3]Vb3 

[A2I  -  [M]V2  +  [CJD^V  +  [K]D2b2  +  l/2ps[C2]V2b2 

f Aj]  -  IC]V2  +  l/2ps[C1]V3b  +  [KjD^V 

[Aq]  -  tK]V2  +  l/2ps[C0JV4 

and  (A-12) 

IB4]  -  (Mc]D2b2 

(B33  -  tMc]DlbV  +  l/2pslC3cJVb3 

lB2J  -  [McJV2  +  l/2ps(C2c]V2b2 

IB  J  -  l/2ps[C  lV3b 
1  lc 

|»J  - 


For  this  study,  three  elastic  modes  were  utilized  (first  wing  bending, 
second  wing  bending  and  first  wing  torsion)  resulting  in  three  general" 
lzed  coordinates  comprising  the  <j(t)  vector.  The  generalized  coordinates 
correspond  to  displacements  of  the  wing  surface  due  to  each  of  the  elast¬ 
ic  modes.  q^(t)  and  q^(t)  are  displacements  measured  perpendicular  to  the 

elastic  axis  in  ft,  and  q  (t)  is  the  rotation  about  the  elastic  axis  in 

3 

radians.  With  jj(t)  being  3x1  and  ^(t)  being  2x1,  the  resulting  [A^] 
and  [B^ J  matrices  defined  in  Eq  A-12  become  3x3  and  3x2  respectively. 

In  addition,  the  aerodynamic  control  surfaces  were  assumed  to  be  rigid, 

massless  plates  resulting  in  [M  ]  being  a  zero  matrix.  Eq  A-ll  can  now 

c 

be  manipulated  so  that  a  state  model  can  easily  be  obtained.  By  equat¬ 
ing  derivatives  with  powers  of  s  (to  revert  back  to  time  domain),  mult¬ 
iplying  through  by  lA^]'1  and  shifting  the  control  terms  to  the  right  hand 
side  of  the  equation,  the  equations  of  motion  become: 

5.(c)  +  la3Jc|('t)  +  [a^jjU)  +  la^j£(t)  +  laQJjj(t) 

"  +  +  tbl^c(°  +  (b0J4c(t) 

where 

[a3J  -  lA4r1lA3] 
la2]  -  (A4J'1[A2] 

i.0j  -  i*4i-‘i*0) 


lb3)  - 

•V 

tbj  - 

'V  ■ 


-tA4l-‘[»2l 

“lA  I”1 [B  ] 
4  1 

“[a  rl  is  j 

4  0 


(A-14) 


Eq  A-13  can  now  be  cast  into  state  space  form  by  making  the  following 


substitutions  (Ref  15): 


X^t)  ■  £(t) 

X2(t)  -  XL(t)  -  ICjJuU) 

X3(t)  -  X2(t)  -  lc2]u(t)  .  (A-15) 

X^t)  -  X3(t)  "  Ic3]u(t) 

X^t)  -  -la()]X1(t)  -  ta1]X2Ct)  -  ta2]X3(t) 

“  la3 ]X^(t)  +  [c4Jji  (t) 

where 

-  Ib3 1 

i^j  -  ib2]  -  i^n^i 

1*3 1  ■  Ibjl  -  la2l tb3 J  -  ta3Jtc2I 

IcAJ  -  [bQJ  -  [ajJl^J  -  la2JIc2J  -  [ ^ J I e3 J 


Employing  these  substitutions,  Eq  A-13  can  be  cast  into  the  following 
state  space  representation  : 


X(t)  -  lA]X(t)  +  lB]u(t) 
Y(t)  -  lCjX(t) 


(A— 16) 


where  the  state  coefficient  matrix  [A]  is  in  companion  form  and  the  con¬ 
trol  coefficient  matrix  [B]  is  a  full  matrix,  that  is: 


42 


IA] 


1 

o 

o 

w 

o 

1 _ 

V 

0  0  10 

c. 

and  [B]  ■ 

2 

0  0  0  1 

C3 

-a  -a .  -a „  -a. 

c. 

0  12  3 

4 

la  this  study,  two  seasors  were  utilized  resulting  la  the  measurable 
outputs  beiog  the  wing  vertical  displacement  h  and  angle  of  twist  or.  The 
resulting  output  vector  J[(t)  is  2x1  and  is  expressed  as  a  combination  of 
the  generalized  coordinates,  that  is  : 


MM 

yjU)  -  ^(t)  (A-17) 

i-1 

where  NM  is  the  number  of  elastic  modes.  The  output  coefficient  matrix 
[C]  then  takes  on  the  form, 


ICl  -  l  C'  :  0  J 


(A-18) 


where 
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-.00189  .03043  .60840 

.70683  .08840  .68922 


and 


0  «*  a  2x9  sub-macrlx  of  zeroes 


which  results  in 


Y(t)  -  [C)X(t)  - 


h(t) 

«(t) 


(A-20) 


. 

E 

l 
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Appendix  B:  State  Model  Formulation/Root  Locus  Program 


Contained  on  the  following  pages  is  a  listing  of  the  program  RLOCDS 
utilized  to  form  the  state  coefficient  matrices  as  a  function  of  free- 
stream  velocity  and  density.  A  sample  input  file  is  also  given  at  the 
end  of  the  Appendix.  The  Pade'  polynomials  used  as  partial  input  to  this 
program  were  formed  using  the  program  PADE  written  by  Thomas  Noll.  To 
run  this  program  on  the  CYBER  750/74,  the  following  commands  should  be 
executed: 


/ATTACH , LINPACK/UN-APPLIB . 
/ATTACH , EISPACK/ UN*APPLIB . 
/LIBRARY ,  UNPACK  .EISPACK. 
/GET .(INPUT  FILE). 

/GET , LGO-RLOCUS /UN-E  710067. 
/LGO, (INPUT  FILE), DATA. 


The  output  is  then  contained  in  local  file  DATA  which  may  be  cataloged 
by  the  user  for  further  use. 
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o  o  oonoooo 


FROG RAM  RL0CUS(INPUT,0UTFUT,TAPE5=1NFUT,TAP£6=0U  TPUT ) 


PROGRAM  USED  FOR  CANTILEVER  OR  PITCH  ANALYSES 
ALL  DIMENSIONS  ARE  IN  SLUGS  ANO  FT 
NMOOE  NUMBER  OF  GENERALIZED  COORDINATES 

NSURF  NUMBER  OF  ACTIVE  CONTROL  SURFACES  < LE  ?> 


32 

41 
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DIMENSION  C ( 4*4 ),BO( 4*4 ) ,B1 (4*4) ,B2 (4*4) .63(4,4) ,AI(  4,4) 
DIMENSION  G(161)  »'»' (161)  »VEL  1(161 >, VGA{ 960 >  , OMEGA (48) 
DIMENSION  AMAT (4d,46),BMAT(48,48  >,BM0(4»4> .8*1 (4 ,4) ,BM2 (4,4  > 
DIMENSION  BM3(4,4>  *B *4(4,4 >,Z4( 4,4) ,Z3( 4,4) ,Z2( 4,4  >,Z1 (4,4) 
DIMENSION  ZO (4,4) ,23 A(4,4) ,Z2AC4 ,4  >  ,Z1A(4,4) ,Z0A (4,4 ) 
DIMENSION  Z3B(4,4)*Z2B(4,4),Z1E(4,4),Z03(4,4),Z1C(4,4) 
DIMENSION  Z0C(4,4),Z0D(4,4) 

DIMENSION  SR2(48),SI2(48),IV2(48  > , F V2 ( 48 > , OMEG ( 4 8) 

REAL  M,K (4,4 ) ,MC(4,4 ) 

COMMON  /PL/  JV,IV,XA(960),YA(960  ) 


COMMON  /VEL/  CM 0(6, fa), CM  1(6,6), CM2 (6 *6), CM 3(6, 6), D MO, OM1, 

1  DM2 , V 

COMMON  /ACT/  AO (4,4) «A1 (4,4 ) «A2(4,4) ,A3(4,4) ,A4(4,4) ,M ( 4 , 4 ) , 
1C0(6,6),C1(6,6),C2(6,6),C3(6,6), IS , IPR 1 , IPR2 ,X SC  ALE 
COMMON  /EIGEN/  A (48, 48 ) , SR (48) , S I( 48 ) , SRI ( 48 ) , SI  1 ( 48  )  , 

1FV 1(48), IV 1(48), SI HZ  (481 
*V  =  720 
IV  =  0 
IS  =.  12 
ISKAX  =  48 
IVMAX  =  15 
ITOTP  =  16*1 VMAX 
1TOTA  =  I$MAX*IVMAX*ITOTP 
CO  33  I=l,ITOTA 
X A ( 1 )  =  YA ( I )  =  VGA(  I)  =  0.0 
CO  41  I=1,ISMAX 
SR ( I )  =  SI(I>  =  0.0 
DO  45  1=1,48 
DO  45  J= 1,48 
A ( I  * J)=Q .0 


AMAT (I ,J )=0. 0 
EMAT (I ,J )=0 . Q 
CONTINUE 


DO  46  1=1,4 
DO  46  J= 1,4 
BMO ( I , J ) =0 • 0 
“Ml ( I,J)=0.0 
cM2 ( I ,J ) =0 • 0 
BM3( I,J)=0.0 
6M4{ I, J) =0,0 
ZO(I,J)sO.O 
Z1(I,J)=0.0 
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Z'2(I,J>  =  0.0 
Z3 ( I  * J)  =  Q*0 
24(1  *J>=0.0 
20  A ( I*J)=0.0 
21 A  ( I  *  J  j  =0  •  0 
22A ( I « J) =0 • 0 
23 A(I*J> =0*0 
46  CONTINUE 

TIME  SCALING  FACTOR 

3ETA  =  0  *0  01 


BET1=8ET  A 
BET2=9ET1*BETA 
BET3=BET2*BETA 
BET4=BET3*BETA 

NPASS  =  1 
URITE(6«3) 

RE AO  IN  OATA  FOR  PASSIVE  FLUTTER  ANALYSIS 

NUMBER  OF  MOOES*  SOLUTION  TYPE  AND  NUMBER  OF  ACTIVE  SURFACES 
PRINT  COMMANDS  IPR1  =  0  WRITES  ALL  INPUT  DATA 
IPR2  =  C  WRITES  ALL  EIGENVALUES 
IPR1  =  1  NO  WRITE 

IPR2  =  1  WRITES  ONLY  SIGNIFICANT  EIGENVALUES 

NSOLTYP  READ  IN  AS  ZERO  SINCE  ACTIVE  CONTROL 
ANALYSIS  IS  DONE  BY  SEPARATE  ANALYSIS 

READ (5*1)  N MODE* NSOLTYP  * NSURF* IPR1 *  1  PR 2 
C  REFERENCE  CHORO*  AIR  OENSITY*  REFERENCE  SEMICHORD  AND  PLOT 
C  SCALING  PARAMETERS  IN  X-CIRECTION  AND  Y-CIRECTION 
READ <5*5 )  BREF*RHO*S*EF*XSCALE*YSCALE 
UR  I TE (6* 14  ) 

IF(NSOLTYP*EQ.O)  WRITE(6*15) 

IF(NSOLTYP*EQ.l)  URITE(6*16> 

IF(NSOLTYP.E(3.2)  WRITE(6*17> 

UR  ITE (6 «  18 )  NM00E*NSURF 
UR ITE(6* 19 )  BREF*RHO*SREF 
00  20  I=l«NMOOE 
CO  21  U=l*NMODE 
MC ( I *J  )  =  0*0 

21  M ( I  *  U )  =  C(I.J>  =  K(I«J)  =  0.0 

20  CONTINUE 

CONST  =  *5*RH0*SREF 
C  GENERALIZED  MASS 

READ  (5*5  )  <M( I«I ) *I  =  l*NHCOE) 

C  GENERALIZED  DAMPING 

READ (5*5  )  (C (I »I )«I=1*NMC0E) 
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C  GENERALIZED  STIFFNESS 

READ (5«5 )  (K (I *  1  )*I=1*NM0DE) 
IF(IFRl.EQ.l)  GO  TO  39 
URITE(6*3> 

WRITC(6*4) 

CO  24  1=1 *NMOOE 

24  WR IT C(£«  2)  ( M ( I • J  > • 0  =  1 «  N  MO  0  E  > 
WRIT£(6«fa) 

00  22  I=1*NM0DE 

22  URITC(6»2>  ( C ( I » J> * J=1 « NFOOE > 


URIT£(6*7) 

DO  23  1= 1*NM0DE 

23  UR  IT E(o»  2)  < KT I ♦ J> *J=1*NMOOE 
39  CONTINUE 

KM  =  NMOOE*NSURF 

C  PAGE  POLYNOMIALS  OF  AEROCYNAFIC 
C  MODES  ANO  ALL  CONTROL  SURFACES 


RE  AO 

(5, 

5) 

( (COII 

«J> 

f  J=1 

*  NM 

REAO 

(5* 

5) 

C  CC1CI 

«J> 

t  J=1 

*  AM 

READ 

C5* 

5  ) 

(  (C2U 

*J) 

t  J=1 

♦  KM 

READ 

<5, 

5) 

( ( C3  ( I 

*0) 

f  J=1 

*  NM 

REAO 

(5* 

5) 

01*02 

00  = 

1. 

0 

IF  ( I 

PR1 

.EG 

.1)  GO 

TO 

81 

URITE(6*3) 

WRITE ( S*  8) 

00  28  1=1  * NMOCE 
CO  29  J=1*NM00E 
WR ITET6  *  9)  I.J 

29  WR ITE (6*10  >  CO C I  * J> » Cl  II *J> * 
28  CONTINUE 
81  CONTINUE 
VELOCITY  VARIATION 


) 


FORCES  FOR  THE  ELASTIC 

1  =  1  *  NM ) 

1=1 »NM> 

1  =  1  *  NM  ) 

1  =  1*  NM  > 


C2( I *J)*C3(I*J>* 00*01*02 


100  CALL  VEL(CO*C1*C2«C3*00*01*02*NM*8REF*CONST«NPASS> 
C ****************************  *  ************************* 
CALL  URITE6<C0*NM*NM> 

CALL  URI TE6<C1*NM*NM) 

CALL  WRITE6<C2*NM,NM> 

CALL  URITE6CC3*NM*NM> 

WR IT  E (6*  87 )  00*01*02 

87  format (//»ix *"00*01*02  ="*3Ei2.5*/> 

IFTV.LT. 0.0)  GO  TO  25 
WR IT  E  <6*  89 )  OMO  *  CM1 , CM2 
89  FORMAT(//*1X*"OMO,OM1*OM2  ="*3E12i5«/> 

CALL  WRI TE6(CM0*NM*NM> 

CALL  URI TEG ( CM  1 *NM  «NM ) 

CALL  WRITE6(CM2*NM«NM> 

CALL  WRITE6<CM3*NK*N") 

CO  30  I=1*NM0CE 
C  FORM  EIGENVALUE  PROBLEM 
CO  31  J=1*NM0DE 
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A4 ! 1 t J )  =  DM2*M!1.J> 

AI!I*J>  =  A4!I,J) 

A3( 1 1 J)  =  0M1*M  < I«J)*0M2*C!I  »J)*CM3!ItJ) 

A2(I»J)  =  0M0*M!I«J)-*DM1*C!I  *J)  +  GM2*K!  I»J)*CN2!  I  »J> 
Al(ItJ)  s  DMO  *C! I »  J ) ♦DM1*K!I *J)*CM1!I«J> 

31  AO ( I » J)  s  OMO*K! I«J)*CM0 i I»J) 

50  CONTINUE 

CALL  URITE4!AOtNKOD£«NMOOE> 

CALL  URl TEA ( A1 »N MODE  «NM OCE ) 

CALL  WRITE*! A2fNM0D£tNM0CE> 

CALL  WRITE4! A3 t NMOQE .NMO CE ) 

CALL  WRITE4! A4,NM00E ♦NMOCE) 

CALL  URI TE4 ! AI *NMOOE  »NMO  CE ) 
IF!NS0LTYP.E6i#2.AND.NPASS.£Q.2>  60  TO  103 
CALL  INVERT(AItNMOOE) 

CALL  WRITE4< AIfNPOnE*NMOCE) 

CALL  MULTMM! AI*A0*80 *NMOCE« NMOOE *NMOCE > 

CALL  MULTMM (AI*A1*B1 *NMOCE*NMODE«NMCCE> 

CALL  MULTMM! AI.A2, 82 ♦NMOCE.NMODE*NMODE> 

CALL  MULTMM  <AI«A3«B3 *NMO CE«N MODE tN MODE) 

CALL  URITE4!BO,NMOOE»NMOOE> 

CALL  WRITE4!B1,NM0CE»NM0CE> 

CALL  URITE4!82#NHGDE.NH0CE) 

CALL  UR1TE4(S3«NM00E«NM0CE> 


CO  50  IsltNMOOE 
00  51  J=1*NSURF 
JU=NM00E*NSURF-J*1 
6M0 ( I«  J) -CMO i l * JU) 

BM1( ItJ)*CMl !I*JJ> 

8M2! I*J)=MC! I«J> *CM0*CM2 (I • JJ> 

6M3 ( I  * J) =MC !I*J)*0M1*CM3!I*JJ) 
EM4(I,J>=MC(I*J)*0M2 
51  CONTINUE 
5  C  CONTINUE 

C 

CALL  URI TE4  !BMO*NMOOEtNSURF) 

CALL  WRITE4!BM1»NM0DE»NSURF) 

CALL  URITE4(DM2*NM0DE«NSLRF) 

CALL  URITE4! SH3*NMQ0E»NSURF> 

CALL  WRI TE4(BM4*NM00E«NSURF) 

C 

CALL  MULTMM !AI»BH4*Z4»NMCDE*NM0DE*NSURF) 
CALL  MUL  TMM  !AI»BM3*Z3*NMCOEtNMODE*  NSUR  F ) 
CALL  MUL  TMM  (  A I  *8iT2  *Z2  »NM  COE  tNM ODE*  NSURF) 
CALL  MULTMM ( A I *0M1 »Z 1,NM00E tNMDQE* NSURF) 
CALL  MULTMM! A  1*0 MO *2 0 »NMOOE *NMOO E*NSURF ) 
C 

CALL  MULTMM! 83 *Z 3*22  A* NM COE* NMODE* NSURF) 
CALL  MULTMM! 62*25*21 A*NMGDE tNMODE* NSURF > 
CALL  MULTHM!B1 *Z3*Z0 A»NMCOE*NMOOE* NSURF) 
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CALL  URI TE4(24*NH0DEfNSURF> 

CALL  WRI TE4{Z3»fcMO0E»NSURF) 

CALL  yRITE4<22»N:«0  3£*NSURF) 

CALL  URI  IE4(21*NM0QEtNSURF) 

CALL  URITE4(20*NMOOE*NSURF> 

C 

CALL  WRITE4<23A*NM00EiNSURF> 

CALL  URITE4(22A»NM0Q£*NSURF> 

CALL  WRITE  4(21A»NMQ0E»NSliRF) 

CALL  WRITE4(20A,NMODE*NSURF) 

C 

C  FORM  8MAT  FOR  STATE  MODEL 

C 

IlsMMOOE 

12=2*NM0CE 

I3s3*NH00E 

C 

DO  42  1= ltNMODE 
DO  42  d=l*NSURF 
Z3B ( I  * J  >  =  -Z3  ( I « J ) 

22B < I » J)  s  -22(1 f J )  «  Z2A(ltJ> 

42  CONTINUE 

CALL  HULTMM(B3*Z2B*Z1C»NM00E*N>«0DE»N$URF> 
CALL  URITE4(Z1C,NM0DE»NSURF) 

00  42  I-1«NM00E 
00  43  jsl,NSURF 

21B(  I* J)  =  -Zl(IfJ)  ♦  Z1 A ( I » J)  -  Z1C(I*J> 

43  CONTINUE 

CALL  MULTMM(B2»Z2B«ZOC»NMODc»NHODE*NSURF> 
CALL  HULTMM (B3«Z16«ZQD*NM0DEtNM9DE*NSURF ) 
CALL  URITE4(Z0C«KMOO£*NSURF> 

CALL  UR1TE4(ZOO«.\MOOE*NSURF) 

DO  44  I=1*NM0DE 
00  44  jsltNSURF 

20B (  I* J)  =  -ZO(I»J)  ♦  20ACI » J)  -  ZOC(I*J> 

44  CONTINUE 
C 

00  53  Is ltNMOOE 
00  54  Js 1 tNSURF 
fiMAT<I,J)s8ETI*Z3B(I»J) 
BMAT<I1M*J>=BET2*Z23< I* u) 

BMAT <I2M» J)=BET3*Z10(  I* J) 
BMAT(I3*I*J)=BET4*Z0B(I» J) 

54  CONTINUE 
53  CONTINUE 
C 
C 

C  FORM  UNAUGMCNTED  DYNAMIC  MATRIX 
Nls  NMO OE *4 
N2  s  NM0CE*3 
.C 
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CO  82  I  =  1  ♦ N 1 
00  82  J= 1 «N1 
A<  I*vi)=0  .0 

CO  36  I=1*N2 
J  =  NMODE*I 
A(I,J)  =  1.0 
K3  =  N2*l 
II  =  0 

00  37  I=N3«N1 

11  =  NMOOE 

12  =  2*NM00E 

13  =  3*NM00E 
II  =  11*1 

JJ  =  0 

CO  38  J=l»NMODE 
JJ  =  JJ*1 

A(ItO)  =  -BOUIiJJ) 

11  =  11*1 

A (  I  *  1 1 )  =  -B1 < 1 1 *JJ) 

12  =  12*1 

At  I* 12)  s  -B2(II*JJ> 

13  =  13*1 

At  I *13)  -  “B 3(11  * JO ) 

CONTINUE 

FORM  AMAT  MATRIX  FOR  STATE  MODEL 

CO  83  1=1, N1 
□0  83  U=1*N1 
AMAT  ( I  * J  )=  0 • Q 

DO  70  1=1, N2 

J=NMOOE* I 

AM AT  1 1  * J  )  =  1.0 


N3=N2*l 
11  =  0 

00  71  I=N3,N1 
I 1=NM00E 
I2=2*NM00E 
13=3  *N.MQOE 
II=I 1*1 
JJ=0 

CO  72  J= 1*NM0DE 
«JJ=J<J*1 

AMAT ( I  * J  >  =  -BO(II*Ow)*BET4 
11=11*1 

AMAT  ( I  « 1 1 )  =  -B 1 < 1 1 «  JJ  > *SET3 
12=12*1 

AMAT (1*12)  =  -B2(II»UJ) *£ET2 
13=13*1 


Ql 


AMATU.I3)  =  -B3<II. JJ) *8ET1 
72  CONTINUE 
71  CONTINUE 

URIT£<6»95> 

95  FORMAT(//tlX*"AMAT  FOR  STATE  MODEL  :"♦/> 

DO  55  1=1 *N1 

55  WRITE(6»98>  ( AMAT ( I » J) » J =1 , 1 0 ) 

WRITE (6*97) 

DO  56  I = 1 « N 1 

56  'JR  I T  E  (  6»  98  )  <  AHA  T  <  I  *J>  •  J  =  1 1  *N1  ) 

C 

UR  IT  E  (  6  *  94 ) 

94  F0RMAT(//*1X»"5MAT  FCR  STATE  MODEL  :■*/> 

00  57  1= 1  * N1 

57  WRITE<6f98)  < BMA T< I* J) « J=1 » NSURF ) 

WR I TE (6  *  97 ) 

C 

CALL  RG(48*Nl»AMAT*SR2*SI2*0»OUMMtIV2»FW2tIERR) 

WR I TE (6*  200 ) 

2C0  FORMAT(//*lXt "EIGENVALUES  OF  AMA T  MATRIX  :«.//»6X, 
&"REAL"*8X,«IMAG".8X* "OMEGA"*/) 

DO  80  I=1«N1 

0MEG(I)=SQRTCSR2(I)**2+SI2(I)**2) 

WRITE(6f 210)  SR2 (I >  *SI2 (  I) tOMEG (I) 

£10  F0RMATC1X.3E12.5) 

80  CONTINUE 
C 

C  PERFORM  EIGENVALUE  SOLUTION  (REQUIRES  EISPACK  ROUTINES) 
WRITE(6,99) 

99  FORM AT  < // *1X  * "  A  MATRIX  (IN  MAIN)  I"t/> 

DO  90  I=1*N1 

90  URITE(6«98)  ( A ( I  * J) . J=1 • 10 ) 

4K ITE (6*  97) 

97  FORMAT!//) 

98  FORMAT<1X*10E12.5) 

DO  91  I=1*N1 

91  URITE(6»98)  ( A( 1  ,  J  )  , J=ll »N1 ) 

CALL  RG(48*Nl*AtSP«SI»0*CUMMy*I  VltFVl* I  ERR  ) 

WRIT£<6. 11) 

PHI2  =  6.28318 
DO  40  1= 1 *N1 

C  CALCULATE  ESTIMATE  OF  STRUCTURAL  0AM PING 
SO  =  0.0 

CMEGA(I)=SQRT(SR(I)**2*SI(I)**2> 

IF(SI(I) .GT.l.E-2)  GO  =  2. *SR ( I ) /OMEGA ( I  ) 

IF<ABS<SI<  I)  ).LT.l.E-2)  SKI)  =  0.0 
SIrt2<I)=SI<I)/PHl2 
JV  =  JVM 
XA(JV)  =  SR (  1 ) 

YACJV)  =  SIU) 

VGA< JV)  =  V 

IF<SI(I> .EG. 3.0. ANO.SRt  I  >.EO.O.O >  GO  TO  40 
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SL. 


1FUPR2.E0.1.AND.SI!  I).E3.0. 

1  60  TO  40 

VRITE!6.2>  SRCI),SI(I)*SIH2(I),G0»0MEGA<I) 

40  CONTINUE 

UR IT£(6* 13) 

NP ASS  =  2 
30  TO  100 

25  CONTINUE 

CO  26  J=1»IT0TA 

IF<TA(J)  .LT.O.O)  XA(J)  *  YA(U)  =  VGAIU)  =  0* .  0 
IF(A3S(YA(J) ) .GT.YSCALE)  XA!J>  =  YA(J)  a  VGA(J) 

26  IF(ABS(XA(J) ) .GT.XSCALE)  XA!J)  =  YA<J)  =  V3A1J) 
IFCNSOLTYP.LE.l)  CALL  PLCT1 
IFINSOLTYP.GE.l)  CALL  PLCT2 
IF(NSOLTYP.EQ.l)  CALL  PL0T3 

1G  =  0 

C  CALCULATE  STRUCTURAL  OAMPING  FOR  PLOTTING 
CO  12  I=1*IT0TA 
YA  (  1 )  =  YA  ( 1 )  *PH  12 
IF(YA(I )  .EQ.C.O)  GO  TO  12 
_  IG  =  IG*-1 

VELl(IG)  a  VGA(I) 

W< IG)  a  SQRT!XA! I)**2+YA!I)**2) 

PS  a  XA ( I)/U< IG) 

GO  a  2.*PS 

OL  =  PHI2*PS/SQRT!1.-PS**2) 

PHI2S  a  PHI  2  **2 

G(IG)  a  2.*PHI2*0L/!PHI2S-0L**2) 
IF(AES(GO).GT.0.5)  G(IG)  a  0.0 
U(IG)  a  U( IG ) /PHI2 
12  CONTINUE 


C  ESTABLISH  ZERO  FOR  FREQUENCY  VERSUS  VELOCITY  PLOT 
IG  a  1G*1 

UCIG)  a  GCIG)  a  VELlilG)  =  0.0 
CALL  PLi)T4(G.U,VELl»IG) 

1  FORMAT (5 15 ) 

2  FORMAT!10X*6E12.5) 

3  FORMAT ( 1  Hi ) 


0.0 

0.0 


4  FORM ATC// * 1 0 X* ‘GENERALIZED  MASS  MATRIX*) 

5  FORMAK6E12.5) 

6  FORMATC//*10X.*GENER ALI2ED  OAMPING  MATRIX*) 

7  FORMAT(//tlOX.*GENER ALIZEO  STIFFNESS  MATRIX*) 

8  FQRHAT(//*1QX*"HPA0E  POLYNOMIALS  FOR  AERODYNAMIC  FORCES  " 
1*(C0*C1*S*C2*S**2*C3*S**3/D0*D1*S*D2*S**2) ") 

9  FORM AT <//. 1 0 X » *P OLYNOMI AL  FOR  COEFFICIENT  C * ♦ 1 1  * *t *♦ II , * ) *) 

1C  FORMAT!10X.4E12.5) 

11  FORMAT !//.10X, "EIGENVALUES  FOR  PASSIVE  SYSTEM  :**,//.12X. 

l"REALCR/S) ".3X,"IMAG(R/S  ) • » 4 X* "I  MAG IHZ > " *4 X » "DAM PING » . 
»4X*"CMEGA(R/S)"»/) 

13  FORMA  I  U0X»3  OH**  ********************  •*...*.*) 

14  FORMAT!/ / *10X**ROOT  lOCUS  FLUTTER  ANALYSIS*) 

13  FORMATUOX.*PASSIVE  ANALYSIS  ONLY*) 


S3 


W  CM 


lc  F0RMAT!1CX»*PASSIVE  AND  ACTIVE  CONTROL  ANALYSES* ) 

17  FORMAT  !10X,*ACTIVE  CONTROL  ANALYSIS  ONLY*) 

18  FORMAT!//, 10X,*NUMB£R  OF  MOOES  =  * ,  12 ,/ 1 3X * *NUMBER  OF* 
1*  CONTROL  SURFACES  =  *»I2> 

19  FORMAT!//. 10 X, *R EFER ENCE  SEMICHORD  =  *,E12*5*  FT*/10X, 
1*AIR  DENSITY  =  *,E12.5,12H  SLUGS /F T**3/l CX *R EFER ENCE  * 
2*SEMISPAN*,£12.5*  FT*=  *) 

STOP 

ENO 

C*********** ******** ********************** **************** 
SUBROUTINE  VEL ! C 0 »C1 »C2 , C3 ,D0 ,01 ,02 ,N,3,C0\ST,NP ASS) 
C********* ****** ******************  ************  ****** ****** 

C  THIS  ROUTINE  ALLOWS  FOR  VELOCITY  VARIATIONS  FOR  THE  RLOCUS. 
C  NEGATIVE  V  STOPS  THE  SOLUTION  PROCESS 

OIMENSION  C0!6*6)*C1!S»6>«C2!6*6),C3!6,6) 

COMMON  /VEL/  CMO  !6  ,6 ) ,CM1 16 ,6) ,C M2 ! 6,6 ) , CM3 ! 6, 6 ) , 

1  CM  0  ,  DM  1 ,  CM  2  ,  V 
V  =  *000  00E  +  00 
IF ! NP ASS .GT . 1 )  REA0!5,1)  V 
IF ! V*LT* 0* )  GO  TO  4 
WRITE! 6,3) 

3  FORMAT ! lhl ) 

UR ITE !6 , 2 )  V 

3  FORMAT!/ »10X, ‘VELOCITY  =  *,£12.5*  FT/SEC*) 

I  FORMAT !E12*5  ) 

00  10  1=1, N 

/§  00  11  J-l ,N 

CMO ! 1 ,J)  =  CO! J,I)*C0NST*V**4 
CM1 ! I *J)  =  Cl!w,I)*C0NST*B*V**3 
CM2 !  I* J)  =  C2!J.I)*C0NST*!B*V)**2 

II  CM3!  I, J)  =  C3!J,I)*C0NST*V*5**3 

10  CONTINUE 

CMO  =  00  *V**2 
0M1  =  01 *B* V 
CM2  =  02*0** 2 

4  CONTINUE 
RETURN 
ENO 

C****» ********************************  ******** 

SUBROUTINE  MULTMM ! A , S,C , N,M,L) 

C* ******************** ********  ********  ******** 

C  THIS  ROUTINE  PERFORMS  MATRIX  MULTIPLICATION  FOR  MATRIX 
C  !  A )  WITH  DIMENSIONS  N,M  TIMES  MATRIX  IB)  WITH  DIMENSIONS 
C  K *L  ANO  STORED  IN  MATRIX  !C)  WITH  DIMENSIONS  N,L 
Cl MENS  ION  A! 4,4) ,8 14  ,4) ,C!4,4) , V (4) 

00  1  1=1, N 
DO  2  J=1»L 
00  3  K=1,M 

TOT  =  T0T*A!I»K)*B!K»J) 

V !  J  >  =  TOT 
00  4  J=1 ,L 
4  C! I , J)  =  VIJ) 
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1  CONTINUE 
RETURN 
ENO 

C* **********************  ********** 

SUBROUTINE  I NVER  T ( A  *  N> 

C* ««««««** ************** ********* 

C  THIS  ROUTINE  PERFORMS  MATRIX  INVERSION  (UNPACK  ROUTINES). 
C  INVERT  MATRIX  (A)  WITH  DIMENSIONS  N •  N  ANC  STORE  IN 
C  MATRIX  (A) 

DIMENSION  A<4*4) *V1(4) , V2(4>  *DET(2> 

CALL  SGEC0(A*4*N*V1*RC0NC*V2) 

CALL  SGEDI(A*4*N*V1.CET*V2»1> 

RETURN 

ENO 

RETURN 

END 


*  **  * 

SUBROUTINE  URI TEA ( A*  N*M> 

i  * 

DIMENSION  A ( 4 *4 ) 

VRITE(6*110) 

CO  1  1  =  1  *  N 

1 

bRITE(6*190)  (A( I*J)«J=1«M) 

100 

F0RMAT(1X*4E12.5) 

110 

FORMATt//) 

RETURN 

END 

♦  *** 

SUBROUTINE  W R I TE6 ( A * N * M ) 

r*  ♦  * 

DIMENSION  A ( 6» 6 > 

WRIT£(6» 110) 

CO  1  1=1 »N 

1  UR  IT  E  (6* 100)  <  A  (  I  *  J  ) t J= 1 *H ) 

100  F0RMAT(IX,6E12.5> 

110  FORMAT (//) 

RETURN 

ENO 

C* ****************************************************  ********** 
SUBROUTINE  PLOT ( X* Y » LX » A « IC *B* HQ  *LL »XQELT* TI TLE » IDM*  TA  PEN ) 


FLOTR*  A  POINT  PLOT  ROUTINE  ADAPTED  FROM  360  PLOT 
FLCTS  A  GRAPH  OF  ONE  OR  MORE  CURVES  FROM  GIVEN  SETS  OF 
RECTANGULAR  COOROINATES 


CALL  PLOTR<X»Y*M,A*IC.B»MP*LL«XDELT*TITLE* IDM) 


NAME  OF  A 
INATES  OF 
INATES  OF 
X( 1*  N)  TO 
THE  CURVE 


2  DIMENSIONAL  AKRA 
ALL  THE  CURVES  TO 
THE  NTH  CURVE*  FOR 
X ( N « N )  WHERE  M  IS 


Y  CONTAINING  THE  X  CCOWO- 
3E  PLOTTED.  THE  X  COORD- 
EXAMPLE*  ARE  STORED  FROM 
THE  NUMBER  OF  POINTS  IN 
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i* 


c 

C  Y 

c 

c 

C  M 
C 
C 
C 

C  A 
C 
C 
C 

c  ic 

c 

c 

w 

C  B 

c 

c 

c 

c 

c 

c 

c 

c 

C  HP 

C 

C 

C  L 

c 

c 

c 

C  XOELT 

C 

C 

C  TITLE 

C 

C 

C  I  OH 
C 

c 

c 

c 

c 

c 

C  TAPEN 

C 

C 


PAS  THE  SAME  OIMENSICNS  AS  X  AND  CONTAINS  THE  Y 
COORDINATES  * 

NAME  OF  ONE  DIMENSIONAL  ARRAY  SET  UP  BY  THE  USER. 

M(N)  IS  THE  NUMBER  OF  POINTS  IN  THE  iM-TH  CURVE*  WHOSE 
FIRST  POINT  IS  AT  X(1,N)  AND  Y  ( 1  •  N  )  • 

NAME  OF  ONE  DIMENSIONAL  ARRAY  SET  UP  BY  THE  USER. 

A(N>  IS  THE  CHARACTER— LEFT  ADJUSTED  Ih  A  A  BYTE  WORD 
TO  BE  USED  IN  PLOTTING  THE  \'-TH  CURVE. 

THE  NUM3ER  OF  CURVES  TO  BE  PLOTTED.  IT  MUST  BE  LESS 
THAN  OR  EQUAL  TO  THE  SECOND  DIMENSION  SPECIFIED  FOR 
THE  X  AND  Y  ARRAYS  IN  THEIR  DIMENSION  STATEMENT. 

1  NO  BORDER*  NO  AXIS 

2  BORDER*  NO  AXIS 

3  NO  BOROER*  X  AXIS  ONLY 

4  BORDER.  X  AXIS  ONLY 

5  NO  BORDER*  Y  AXIS  ONLY 

6  BORDER*  Y  AXIS  ONLY 

7  NO  BOROER  «  BOTH  AXES 

8  BORDER*  BOTH  AXES 

0  SINGLE  PAGE  PLOT  DESIRED 

1  MULTIPLE  PAGE  CR  FRACTION  OF  PAGE  PLOT  DESIRED 

1  PLOT  GIVEN  POINTS 

2  SEMILOG  (PLOT  LOGS  OF  Y  COORDINATES) 

3  LOG-LOG  (PLOT  LOGS  OF  X  ANO  Y  COORDINATES) 

0  INDICATES  DELTA  X  IS  TO  BE  CALCULATED 

OTHERWISE  SPECIFY  OELTA  X  IN  FLOATING  POINT. 

NAME  OF  THE  ARRAY  IN  WHICH  THE  TITLE  TO  HEAD  EACH 
PAGE  IS  STORED.  TEN  WORDS  ARE  ALWAYS  PRINTED. 

THE  FIRST  DIMENSION  SPECIFIED  FOR  THE  X  ANO  Y  ARRAYS 
IN  THEIR  OIMENSI CN  STATEMENT.  IT  MUST  3E  AT  LEAST 
EQUAL  TO  THE  NUMBER  OF  POINTS  IN  THE  CURVE  WHICH  HAS 
THE  GREATEST  NUMBER  OF  POINTS  OF  THE  CURVES  TO  BE 
PLOTTED. 


OUTPUT  TAPE  FOR  PLOTR 


DIMENSION  X(1)*YU)*LX(1>»A(1>*PLINE(11>.XAX(6)*TITLE<12> 
INTEGER  A»SS  «PLINE  »e .BLNK.PLUS  *T  APEN 
CATA  PLUS»MINUS* II *BLNK  / 1H ♦ * 1H *  * 1 H  I  * 1  OH 
11  FORMAT ( 1 H*  »  5X.12A10  ) 
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12  FORMAT <  1  7H  SCALE/  INCH  X  =  ♦  1PE1  4.  7  * 3H  Y  =  ,  1PE1  4  •  7 .  1 0  X . 

."♦OR-  TOLERANCC/POINT  X  =  "  «  1PE1 A  .  7  . =5H  Y=.1PE14*7) 

13  FORMATd  1X.11A10  > 

14  FORMATd  X.E9.2.1X.11A10) 

16  FORMATdlX.5  <1H  +  .19X),1H'*/7X,5CE9.2.11X>  *E9.2) 

16  FORMAT d 7H  MINIMUM  X  =  t 1PE1 4. 7 .5H  Y= ♦ 1PE 1 4 . 7 . 23 X ♦ 

."MAXIMUM  X  =  " * 1PE14 .7*  5  h  Y=.1PE14,7) 

XMAX  =  Xd  ) 

YMAX=Y<1 ) 

X0ELT=ABS(X3ELT> 

MP=MQ 

IF (MP.EQ.l.AND.XCELT .EQ.C.JMPsO 
C  INITIALIZATION  IS  NOW  COMPLETE 

C  NOW  CURVES  WILL  BE  SEARCHED  FOR  XMA X » X M IN . YM A X  * YM IN  * AND 
C  XOELT  =M INIMUM  DISTANCE  BETWEEN  ANY  TWO  POINTS  ON 
C  ANY  CURVE 

GO  T0(3.2.1>,LL 

1  XMAX-ABS (XM A  X ) 

2  YMAX=A5S< YM AX ) 

3  XMIN=XMAX 
YMIN=YHAX 
L=2 

IFCXCELT  .NE.  0.)  GO  TO  5 

4  L=1 

5  00  400  J=1.IC 
N=1*I0M*  <J-1 > 

XCHAX=X<N> 

YCMAXsY(N) 

GO  T0(22. 21.20). LL 

20  XCMAXSABS(XCMAX) 

21  YCMAX=ABS(YCMAX) 

22  XCMINsXCMAX 
YCMIN=YC«AX 

GO  TC (23  .27 ) «L 

23  XCOELTsO 

27  NEL=LX(J) 

CO  300  1=1. NEL 
N=I*IOM* (J-l ) 

XTEMP=X(NI 

YTEMPsY(N) 

GO  TC(30 .29.28) .LL 

28  XTEMPsAdS(XTEMP) 

29  YTEMP=ABS( YTEMP) 

30  XCMAX=AMAX1 ( XTEMP . XCMAX > 

XCMINSAMIN1 (XTEHP.XCMIN) 

YCMAXSAMAXK  YTEMP.  YCM AX  ) 

YCMIN  =  AM  INKYTEMP.YCMIN  > 

GO  TC(31  .300  >  «L 

31  IF ( I  .GE.  NEL >  GO  TO  300 

311  M=I«-1 

00  100  K=M.NEL 
N=K+  IOM* (J-l ) 
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CO  TC(32  *32*33) *LL 

32  XDIFF=XTEMP-X<N> 

GO  TO  335 

33  XDIFF=XTEMP/ABS< XC N> > 

IF (XD1FF  .EQ  •  0.>  GO  TO  100 
332  XOIFF=ALOG10 <XOIFF) 

335  XOIFF=ABS(XDIFF) 

IF<X0IFF  *£Q •  0.)  GO  TO  100 

336  IF  IXCOELT  .NE.  0.)  GC  TO  338 

337  XC0ELT=X3IFF 
GO  TO  100 

338  XCOELT=ArtINl (X01FF ,XCOELT) 

100  CONTINUE 

300  CONTINUE 

XMAX=AMAX1(XMAX*XCMAX> 

XMIN=AMIN1CXMIN*XCMIN> 

YMAX=AMAX1 ( YMAX*YCHAX) 

YMIN-AMINItYMIN»YCMIN) 

GO  TO<34,400)*L 

34  IF ( X DELT  .NE.  0.)  GO  TO  346 

345  XOELT=XCOELT 
GO  TO  400 

346  XDELT=AMIN1( XOELT, XCCELT) 

400  CONTINUE 

C  IF  XOELT  IS  STILL  ZERO  THEN  THE  CURVES  ARE  PARALLEL  TO  THE 
C  .  Y  AXIS  SO  XMAX.XMIN  AND  XOELT  ARE  CALCULATED  3Y  OTHER  KEANS 
(•  C  DESIGNED  TO  CENTER  CURVES  CN  A  SINGLE  PAGE 

IFCXCELT  .NE.  0.)  30  TO  365 

347  IF(XPAX) 354*353*352 
*53  XHAX=10.*XMAX 

XMIN=XMIN-10 . 

GO  TO  348 

352  XMIN=XMIN-XMAX 
XMAX-2.*  XMAX 
GO  TO  34 £ 

354  XHAX=XMAX-XMIN 
XH  IN  =  2  •  *  XM  IN 
GO  TO  348 

365  IFCXKAX  .NE.  XMIN)  GC  TO  3365 
3366  XMIN  =  X!TIN-54.*XDELT 
XM AX =XMAX+54.* XOELT 

348  f*Psc 

C  IF  THE  CURVES  ARE  PARALLEL  TO  THE  Y  AX]S*A  SIMILAR  PROCEDURE 
C  IS  FOLLOWED  FOR  YM4X, YM IN* ANOYDELT 
3365  IFIYKAX  .NE.  YMIN)  GO  TO  366 

349  IF(YMAX) 36 2*563*364 

363  YMAX  =  10.  -*YMAX 
YMINsYMIN-10 . 

GO  TO  351 

364  YMAX=2.*YMAX 
YM IN  =  YM I N-YM  AX 
GO  TO  351 
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262 

YMIN=2.*YMIN 

YMAXsYMAX-YMIN 

251 

MP  =  0 

9 

266 

GO  T0(37,358 ,35) , 

LL 

35 

IF  <XMAX  .EQ.  0.) 

GO 

10 

356 

255 

XMAX=ALOG10(XMAX> 

256 

IF ( XrtIN  -EQ.  0.) 

GO 

TO 

338 

25  7 

XMIN=ALOG10(XMIN> 

258 

IF  C Y  MAX. EG • 0  •  )  GO 

ro 

359 

36 

YMAX=ALOG10  C  YM AX  ) 

• 

259 

IF ( Y M IN. EG. 0  •  )  GO 

TO 

37 

361 

YMIN=ALOG10( YMIN) 

37 

YRA,NGE=YMAX-YMIN 

xrange=xmax-xmin 

PRANGE=XDELT*103. 

IF 

the  single  page  option 

IS  DESIRED  BUT  XDELT  IS  TOO  LARGE 

C  FOR  THE  RANGE  OF  X  VALUES  TO  FIT  ON  ONE  PAGERS*  A  NEW  XCELT 
C  MUST  9E  FOUNO 

IF  (MF  .NE.  0.)  GO  TO  376 
274  IF (XRANGE  .EQ.  PRANGE)  GC  TO  376 
375  XDELTsXRANGE/108. 

PRANGEsXRANGE 
276  NPAGE=XRANGE/PRANGE+I. 

Y0ELT=YRA.NGE/50. 

YSF=YRANGE/8.33 

XSF=PRANGE/10.3 

XTP=XDELT/2. 

YTP  =  YOEL  T/2 • 

PXHAX=XMIN*XTP 

C  NOW  THE  PLOT  IS  FORMED  A  LINE  AT  A  TIME,  SEARCHING  EACH 
C  CURVE  ONCE  FOR  EVERY  LINE  CN  EVERY  PAGE  AND  PRINTING  OUT  EACH 

C  LINE  AS  SOON  AS  IT  IS  FORMED 
35  DO  900  K=1.NPAGE 

IF (PXHAX-XMAX  .GE.  XTP)  GO  TO  73 

40  IFIPXMAX  .LE.  XHAX.  GO  TO  401 

404  IF(K-1)3 75 ,375*73 

401  PXMIN=PXMAX-XOELT 
PXMAX-PXMAX+PRANGE 

IF ( K  .NE.  NPAGE)  GO  TO  402 
403  PXMAX=XMAX*XTP 

402  GO  T0<42 *42,41 )»LL 

41  PMXMNS10 .**PXMIN 
PNXMX=1Q.**PXMAX 

42  CONTINUE 

UR ITE (6 , 6000  ) 

6000  FORMAT ( 1H1 , / / ) 

URITE(6, 11 >  (TITLEII >*1=1,12) 

URITE(6,6010) 

6010  FORMAT ( 1  HO ) 

WRITE! TAPEN, 16)  XM  IN  ,YM  I  N, XM AX  * Y MAX 
WRITE  CTAPEN* 12)  XSF , YSF , XTP  *  YTP 
GO  T0(44, 43, 43, 43, 44, 43, 43, 43), 5 
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4!  XTEMP=PXMAX-XTP 

JJJ=(AMIM(XTEHP»XMAX>-PXM1N)/XDELT 
44  YU8= YMAX ♦YTP 
YLB= YMAX-YTP 
CO  700  LLINE=l*51 
LINE=52-LL1NE 
00  450  1=1*11 


450 

PLIN£( I ) =8LNK 

GO  T0(46»45*45)*LL 

45 

YNUB  =  10.  <**YUB 

YNL8=1C. ** YLB 

46 

GO  TC(51. 47*511*47*519* 

47*511*47). fl 

47 

fc=JJJ/10 

NCHAR=MOO(JJJ*10) 

IF (LINE- 1)50  *49*48 

4e 

IF (LINE  .NC.  51)  GO  TO 

50 

49 

IF(N  .LE.  0)  GO  TO  549 

5549 

CO  S00  I =1  * N 

500 

PLIME(  I)  =10  )(♦♦♦♦♦♦♦♦♦♦ 

549  NCHAR=NCHAR-H 

CO  550  J  J=1 « 6CH  AR 
d=JJ-l 

550  CALL  PACK(PLIN£(N*1) *PLUS*J) 

60  TO  5118 

50  CALL  PACK(PLINE(1) «PLUS*0) 

CALL  PACK(PLINE(N+1) *PLUS*HCHAR> 

5118  GO  TC<51. 51*511*511* 519*519*511*511)  *8 

511  IF(YUB)512, 513*514 

514  IF ( YLB  •  GE.  0.)  GO  TO  512 
513  N=JJ J/10 

NCHAR=MOO(JJ«J*10  )  ♦  1 
IF <N  .LE.  0)  GO  TO  5515 
5513  CO  515  1=1  »N 
515  PLINE(I)=10H********** 

5515  DO  521  I=l*NCHAR 
d=I-l 

521  CALL  PACK(PLINE(N+1)  *HI.MUS*J) 

512  GO  TO  (51*51*51*51*519*519*519*519)  ,3 
519  IF (PXMAX )51«516*517 

517  IF (PXMIN  .GT.  0.)  GO  TO  51 
516  I=-PXMIN/XOELT 
NCHAR=MO 0(1*10) 

I=l/10*l 

CALL  PACK(PLINE(  I)*H*NCHAR) 

51  CO  600  d=l*IC 
NEL=LX(d> 

CO  600  I =1 «NEL 
R=I*10M* (J-l ) 

GO  TO (54  *52  *  52 ) t  LL 

52  YTEKP=A8S(Y(M> ) 

IF(YTEMP-YNLB)60  0*58  *53 
IF(YTEMP-YNUB)5<3*600 *600 
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54  IF(Y(M)-YL5>600«56.55 

55  IF(Y(M).GE.YU8>  GO  TO  600 

56  XTEf1P=X(M) 

IF (X  TEMP -PXM IN >6 00 *6  2 *57 

57  IF(XTEMP~PXMAX)62*62*600 

58  GO  T0(56*56*59>*LL 

59  XTEHP=ABS(X(M)> 
1F(XTEMP-PNXMN)600*61*60 

6C  IF(XTEMP.GT.PNXKX)  GO  TO  600 

61  IFCXTEMP.EO.O.)  GO  TO  62 
611  XTE.MP=ALOG10  (XTEMP) 

62  JJs(XTEMP-PXMlN) /XOELT 
GO  TC(56t64,63) «LL 

62  IF(X(M).LT.Q.)  GO  TO  65 

64  IF ( Y  (H)*GE»0»)  GO  TO  66 

65  SS=1 HN 
GO  TO  67 

66  SS=A(J) 

67  NCHARsMOD(JJ.lO) 

N=JJ/10*1 

CALL  PACK(PLINE(N)*SS*NCHAR> 

600  CONTINUE 

66  IF  (M00(LINE*6)  .NE.l)  GO  TO  69 
9C1  IF  (LINE  .NE.  1 >GC  TO  71 
902  WRITE(TAPENtl4>  YMIN  *PLI NE 
GO  TO  715 

65  IF (LI NE  .EQ.  51)  GO  TO  71 

70  WRITE(TAPEN*13)  PLINE 
GO  TO  715 

71  YTEHPsYLQ+YTP 
URITE(TAPEN*14)  YTEHPtPLINE 

715  YUB=YU8- YOELT 
YLB=YL8-Y0ELT 
700  CONTINUE 

72  CO  600  I~l«6 
J=I-1 

XAX( I)=PXMIN+XTP>XDELT*FLOAT(J)*20. 
6  CO  CONTINUE 

WRITE  CTAPENt 15)  XAX 
900  CONTINUE 
72  RETURN 
END 

SUBROUTINE  PACK (WORD »CHAR*IP0S> 
DIPENSION  SPLIT ( 10 ) 

OATA  SPL IT/1  0 *1 H  / 

CEC0CE(10«2*W0RD>  (SPLIT  (I ) «Isl  ,10  ) 
SPLIT(IPOS+l)=ChAR 

£NCOOE( 10*2* WORD)  (SPLIT  (I) «  I=l*10) 
RETURN 

2  FORM  AT ( 1 OA 1 ) 

END 

SUBROUTINE  PLOT1 
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c 3. 


Cl  MENS  ION  S( 1), TITLE (12 ) *X ( 2  40 ) • Y( 24  0) 

COMMON  /PL/  JV*IV,XA(960),YA(960 > 

CATA  S  / 1H*/ * T I TLE  /120HROOT  LOCUS  OF  PASSIVE  SYSTEM 

1 

2  / 

J  =  0 

CO  1  1=721*960 
J  =  d»l 
X(J)  s  X  A ( I ) 

1  YCJ)  =  Y  A( I ) 

CALL  PLOT(X*Y,240*S, 1*7, 0,1,0, TITLE, 24 C *6 > 

RETURN 

ENO 

SUBROUTINE  PL0T2 
DIMENSION  S ( 1 >  *  T  I  TLE  (12  > 

COMMON  /PL/  JV»IV*XA(960),YA(960) 

CATA  S  /1H*/*TITLE  /120HROOT  LOCUS  OF  AUGMENTED  SYSTEM 

1 

2  / 

CALL  PLOT(XA*YA*720«S*1* 7* 0*1*0, TITLE* 72 0*6) 

RETURN 

END 

SUBROUTINE  PL0T3 
DIMENSION  S( 1) *T ITLE (12) 

COMMON  /PL/  JV*IV,XA(960 )»YA(960  ) 

DATA  S  /1H*/, TITLE  / 12  0H COMB INED  ROOT  LOCUS 

1 

2  / 

CALL  PLOT(XA,YA* JV,S *1,7 *0*1, 0*T ITLE, JV,6) 

RETURN 

ENO 

SUBROUTINE  PL0T4(G*W*V*N> 

DIMENSION  G( 161) *W(161) *V(240),S (1) , TITLE 1 (12> .TITLE 2(12) 
DATA  S  / 1H*/ *TI TLE1  /120  FSTRUCTUR AL  CAMPING  VERSUS 

1  VELOCITY 

2  / 

CATA  TITLE2  /120HFREUUEN CY  VERSUS  VELOCITY 

1 

2  / 

CALL  PLOT(V»G,N*S,1,7*0*1»0,TITLE1*N*6) 

CALL  PLOT(V,y,N*S, 1*7*0, 1* 0  * TITL E2 , N*6> 

RETURN 

ENO 
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3  C  2  0  0 

1.  .002378  1.  1000.  1000. 

.326832  .C25608  .12 

0.3 

6.062*  161.23  2649.8  • 

-.285S9E*r0  .77822E-01  -.100S4E*01  -.41ft21F.-Cl  -.«9“17E-02  -.12356E-01 

-.819C0E-'2  -.312 57E  «0 1  -.15436E.Q0  -.664S4E-01  -.13K46E*02  .33?‘>1E«C1 

-.46«69E*',2  *.20996E*01  -.32022C*0Q  -.20622E*C0  -.30?:  16E*C0  — . ? f- 7 07 r» Ci  1 
-.11302E*iri  .25951E-01  .34344f>01  -.16266C*01  .43709E-C1  .2J344E*3n 

.94072E*C0 

•10 3 J5E*C1  .94280E-01  .30429E-C1  .16542E*C0  .13270E-Q1  -.16317E*01 

•  1 0079E*  r  1  -.498 20E  »01  -.171668*00  -.143721*00  -.15f2PE*92  .5f048E*91 

-.39928E* r2  -.15323E*01  *.90763E*00  -.353148*90  -.27644E*on  -,23797E*91 

-.11 7E9E*  *1  .17966E-01  .36789£*01  -.17592E*G1  .21653E-01  .1  71  4C-E*0  0 

• IS 132E-T1 

. 12C65E*  01  .33364E-01  .2965CE*01  .1S8C5E*00  .45343E-01  -.606°fcE*90 

. 12636E*r 1  -.22004E«01  -.38022E-01  -.969918-01  -.36239E*01  .3?487E*01 

•  12790E*? 1  .27946E  *0  0  -.66577E.00  -.17243E*00  .S5C24E-01  -.35^75E*09 

•.58786C-<“2  -.92806E-02  .75659E*00  ».42457E*C0  -.55406E-01  .128«5E-01 

•  60leiE*<,0 

.28137E*fC  -.11552E *00  .34789E*00  .65191E-02  .31736E-01  -.13467C*00 

•  32762E* "9  -.36CS2E *0 0  -.31621E-02  -.20013E-C1  .lb02?E*00  -.1S2S5E+00 

.87P69E*nC  .70012E-02  -.13154E-01  .77903E-02  -.2679PE-02  .134C1E-01 

•12624E-C2  .70142E-03  .54399E-01  -.262S0E-G1  .31042E-01  .24164E-02 

. 34408E-01 

•  10 4 18E* 01  .2S446E  02 

94. 

115. 

1S6. 

187. 

203. 

219. 

-1. 
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Contained  on  the  following  pages  la  a  listing  of  the  program  STMOD 

which  conducts  several  state  model  manipulations  at  the  users'  discre- 

« 

tlon.  The  input  matrices  vary  depending  on  the  option  selected.  Cur¬ 


rently,  there  are  five  selectable  options  available  to  the  user  which 
are  as  follows: 


1)  Singular  Value  Decomposition  and  LU  Decomposition 
of  the  Observability  Matrix 

2)  Singular  Value  Decomposition  and  LU  Decomposition 
of  Che  Controllability  Matrix 

3)  State  Transformation  of  a  Matrix  to  Diagonal  Form 
(Eigenvalues  and  Eigenvectors  Calculated) 

4)  Enter  the  Feedback  Cain  Matrix  C  and  Form  A  +  BG 
Matrix  and  Find  Closed  Loop  Eigenvalues 

5)  Do  Same  as  4  but  In  Addition  Read  in  Observer  Gain 
Matrix  and  Find  Observer  Poles 


The  state  coefficient  matrices  can  be  input  manually  (and  are  prompted 
for  by  the  program)  or  read  from  a  data  file  attached  as  TAPE8.  To  ex¬ 
ecute  this  program  on  the  CYBER  750/74,  the  following  commands  are 
necessary: 


/ATTACH , IMSL4/UN-APPLIB. 

/LIBRARY, IMSL4. 

/GET ,TAPE8-(DATA  FILE  NAME-NOT  NECESSARY  FOR  MANUAL  INPUT) 
/GET, LG0-STM0D/UN-E7 10067. 

/LGO. 


oooooooononoooo 


PROGRAM  STM0L;(  lNPUT  =  /80  »  OUTPUT  =  /feO  .  I)  AT  A= /8  0  *  T  APE  5=  IN  PU  T* 
&TAPE6=0UTPUT*TAPE7  =  DATA*  TAP£8  =  /8  0> 


THIS  PROGRAM  IS  USED  TO  PERFORM  SEVERAL  MANIPULATIONS  OF  A 
SYSTEM  IN  STATE  VARIABLE  FORM  (USING  IMSLA  LIBRARY  RCUTJNFS) 

I.E.  XDOT  =  AX  ♦  BU 

Y  =  CX 


WHERE 

A=STATE  COEFFICIENT  MATRIX 
B=CONTROL  COEFFICIENT  MATRIX 
C=CUTPUT  COEFFICIENT  MATRIX 


REAL  A(12*12>*B(12.12)*C(12*12)*CA(12*12).GC12*12)»G1C12*12) 
REAL  P(12«12)tUTl(12«12> *V< 12*12 )* VI (12,12 )*CAl( 12 ♦ 12) 

REAL  VIK(24)*UT(12«12)*S(12)tBA(12*12)*3Al(12*12)«Pl(12tl2) 
REAL  LUC 12*12)*EGUIL(12)«LC12*12)*U(12*12)*'JK3(1H0) 

REAL  AT(12*12)*'jlKl<lb8)  *TREAL(12»12)*TINVR<12*12  ) 

REAL  ERC12)*EI(12) *ETR< 12) *ETI( 12) * TRIO (12*121 *TLU <12*12 ) 
REAL  EVMAG<12>  »WK2< 12), ATP< 12*12  )*TINVBR C 1 2* 12 > « TTRA N ( 12 * 12 ) 
REAL  G<1 2. 12 ) *9G <12. 12 >*ABG( 12*12) .UKAC 12) *BGT I C 12*12) 

REAL  T9X9<12.12>  .TT9 X9 < 1 2. 1 2 > * TTT9 ( 1 2, 12  )  *  AM  AT ( 12* 12 > 

REAL  CTREAL(12*12),K1<12*12) *KCT (1 2  *  12 ) * AE ( 1 2* 12 ) *  AC  1 ( 1 2 . 12 ) 
REAL  AE2(12*12),AE2A (12*12) .2REALC12.12) .CTTC12.12) 

REAL  KIT (12*12) *89 < 1 2* 12 ) . A9 ( 12* 12) ,T9  ( 1 2  *  12  >  *  TT  T  ( 12  *1  2  ) 

REAL  T9 I C12*12)»6G9<12*12) *ABG9< 12*12) .A91C12.12 ) 

REAL  A9CL(12*12) *G9 < 12* 1 2) * BMAT ( 12  *  1 2 ) 

INTEGER  IPVTC12) 

COMPLEX  EVALC12) *EVALTC 12) *2 (12*12 ) *TT( 12*12 )*2N*ZNT 
COMPLEX  AC(12.12)»T1 ( 12  *  12 ) * T< 12 *1 2 ) .ZTTC12.12) 

COMPLEX  TINVC<12*12)  ,WA ( 168)*TID (12*12) *ZTEHPl 12 *12) 

COMPLEX  BC( 12*12  )*TINVB( 12*12) »BZ( 12*12) 

COMMON/ I N OUT /NI N*N OUT *M OUT 
C 

C  INITIALIZE  ARRAYS 

C 

CO  1  1=1*12 
1PVT  (  I )  =  G. 0 
EQUI L ( I ) -0  «  0 
S( I ) =0.0 
ERCl )=0 . 0 
E I ( I  )  =  0  •  0 
ETR  <  I  )  =  G  •  0 
ETI < I >=0 .0 
EVMAGC I ) =0 .0 
WK2C I )=0 .0 
VK A  (  1  >  =0  •  0 

F2(l)=CMPLX(0. 0*0.0) 


CJL 


EVALCI)=CMPLX<0*0*0»0) 

EVALT<I)=CMPLX<0.0*0.0> 

CO  1  J=l,12 

A  ( I »  o )  =0  .0 

B( I* J)=0  «0 

C(I»J)=0.0 

CA(I » J )  =  0 • 0 

G(I»J)=0.0 

P(ItJ)=0.0 

UT1( I»J)=0.0 

V(hJ)sO.O 

LT(I*JI=0.0 

S1(I»J)=0.0 

B  A  ( I  »  J  )  =  0  •  0 

BAH  I  ♦  J)  =0 . 0 

PI (I *J)=0.0 

LU(I *J)=0.0 

L(I»O)=0.0 

U<lt«’>=0.0 

TLU ( I  » J) =0 • 0 

TREALCI. J)=0.0 

TINVRU*  J)  =  0.0 

TRIO  <  If  J  )s0 • 0 

ATR(  I » J) =0 . 0 

TINVBRCI * J )  =  0 • 0 

TTR AN ( I  *  J) =0 • 0 

G<I» J>=0 .0 

EG ( I »J>=0.0 

ABG ( I  * J  >  =0 . 0 

BGTI < I»J)=0. 0 

TTT( I*J) =0 • 0 

T9X9(I* JJsO.O 

TT9X9 (I  * J)  =  0 *0 

TTT9(I»J)=0.0 

AHAT ( I » J )=0 . 0 

EHAT  ( I )  =  0 • 0 

CTREAL<I*J)=0.0 

K1(I»J)=0.0 

KCT<  I* J)=0.0 

AE(I *J)=0.0 

AE1C  If J)=0.0 

AE2  ( 1 1 J ) -0  *0 

AE2A(  If J)  =  0.0 

ZREALdf  JJsO.O 

CTT ( If J>=0.0 

K1T( If J) =0.0 

P9<IfJ)=0.0 

G9(If J)=0.0 

A9  ( I «J)  =  0.0 

T9(If J)=0.0 

T9K  If  J)  =0  •  3 

£G9< If J)=0.0 
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ABG9 ( I  * J )  =  0 . 0 

A91  ( I  *  J) =0*0 

A9CL ( I . J  )  =  0 • 0 

BC(I  *J)=CMPLX(Q.O.  0.0 

TINVDU* J)=CMPLX  (0.0 .3.0  > 

ZT£MP(I*J)=CHPLX(Q.Q«O.G) 

TINVCCI. J>=CHPLX(0.C.0.0> 

TIO( I* J) =CMPLX(0.0 tO.O) 

Z(ItJ)=C.*tPLX<0. 0,0.0) 

TT(I.J)=CMPLX(0. 0.0.0) 

ZTT ( I»J)=CMPLX(0.0»0.0) 

AC(1 1 J)=CMPLX(0.0.  0. 0) 

T1 <1 ,J)=CMPLX( 0.0*0. 0) 

T(I,J)  =  CMPLX (0.0*0. 0  ) 

1  CONTINUE 

DO  2  1  =  1.29 

2  WK ( I ) =0 • 0 

C  SET  UT  EQUAL  TO  IDENTITY  MATRIX  ON  INPUT 
CO  3  J=1 .12 

3  liT  (J»J)  =  1.0 
CO  4  1=1 .168 
WX1 ( I )=0 .0 

WA(I )=CMPLX(0. 0.0.0) 

4  CONTINUE 
NIN=5 
N0UT=6 
POUT  =7 
KIN=8 

JP ASS= 1 

WR I TE (MOLT .3 10 ) 

310  FORMAT ( 1  HI  I 

URITE(NOUT.IOC) 

100  FOR .4 AT(/. IX. "ENTER  N .  THE  OROER  CF  THE  SYSTE*  >■  > 
REAOtNIN.*)  NSYS 
WRITE(NOUT.llQ) 

110  FORMAK/. IX. "ENTER  M.  THE  ORDER  CF  THE  CONTROL  >") 

READ (NIN  **  )  MC 
C* ..*****. *•*..*••*****•*«* 

N=NSYS 
*  =  MC 

C* ******** ••*»*•** *•«*«.*•• 

UR ITE ( NOUT* 1 02 ) 

102  FORMATC//.1X."*****  MENU  -  ENTER  ONE  CHOICE  :  *****",//. 1 X. 

S"<1)  SINGULAR  VALUE  DECOMPOSITION  ANC  LU  DECOMPG  ?  I  T  I  ON  "  .  /  * 
SIX."  CF  OBSERVABILITY  MATRIX  G".//. 

SIX. "(2)  SINGULAR  VALUE  CECOMPOSI TION  A NO  LU  DECO vPOS  IT  ION". 
S/.IX."  OF  CONTROLLI3 I  LI TY  MATRIX  P"»//.l<, 

S " ( 3 )  STATE  TRANSFORMATION  OF  A  MATRIX  TO  DIAGONAL  FOcM"»/ 
S.1X,"  X  =  TZ  (EIGENVALUES  AND  EIGENVECTORS  CALCULATED)" 
A. //.IX."  (4)  ENTER  GAIN  MATRIX  G  AND  FO*M  A.*3  MATRIX  ANC" 
S./.1X."  FIND  CLOSED  LCOP  EIGENVALUES".//* 

>.lx*"(5)  DO  SAME  AS  <**)  jLT  IN  ADDITION  READ  IN"./. 
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WlX*"  OBSERVED  GAIA!  MATRIX  ANO  FIND  OBSERVER"* 

& "  EIGENVALUES"*//) 

READ (N1N  «* )  ICHOICE 
WR IT  £ (NOUT, 1 13 ) 

113  F0RMAT(//*1X, "CO  YOU  WANT  TO  ENTER  STATE  MATRICES  (ENTER  1>" 

1  */  » 1  X  ,  "0  R  DO  YOU  WANI  TO  READ  THEM  FROM  A  FILE  (ENTER  2) 

&  ?"♦/> 

RE AD (NI ft**)  KREAO 
IF (KREAQ  .NE  *  2 )  GC  TO  9 
REWIND  KIN 

CALL  RE AOF ! A  »N  *W  »K IN ) 

WRIT E (MO UT* 2 05) 

CALL  WRI rE(A*N,N*MOUT) 

IF(ICHOICE.EQ.l)  GO  TO  8 
CALL  READF(S,N*M«KIN) 

WRITE! MO UT*210) 

CALL  URITE(6*N,M*MQUT) 

IF(ICHOICE.LT.R)  GO  TO  1* 

CALL  READF(G*M«N*K1N) 

WRITE(M0UT*212) 

212  F0RMAT(//«1X, "FEEDBACK  GAIN  MATRIX  G  I",/) 

CALL  WRI TE(G,M,N*MOUT> 

IFUCHOICE.LT. 5)  GO  TO  13 

CALL  REA0F(K1,N,M*KIN) 

WRITE (MO UT« 2 13) 

213  FORMAT!//, IX, "OBSERVER  GAIN  MATRIX  K  :",/> 

CALL  WRIT£(K1*N,M,M0UT) 

GO  TO  13 

8  CONTINUE 

CALL  REABF(C,M,N*KIN) 

WRITE(MOUT*215> 

CALL  URITE(C,M,N,HOUT) 

30  TO  13 

9  CONTINUE 
WRITE!  NO UT *115) 

115  FORM  AT!//, IX , "USE  LIST  DIRECTED  READ  FOR  THE  ARK  AYS"  ,/ « IX, 

& "ENTER  I*J*NON-ZERO  ENTRY  A(I*U)  ",//> 

10  WRITE! NO UT*120) 

120  FORMAT!/ ,1X»"ENTER  A  MATRIX  !N  X  N)  S",/) 

CALL  REAO(A,N«N,NIN) 

CALL  W  R I TE ( A , ft  *  ft  *ft  OU  T ) 

WRITE! NOUT  *200) 

200  FORrtATUX, "ENTER  0  TO  ACCEPT*  1  TO  CHANGE  >") 

REAO  !NIN  »* )  IANS 

IF(1 ANS.EO.l )  GO  TO  ID 

WRITE(MOUT,205> 

205  FORM  AT!/ *1X* "A  MATRIX  :",/) 

CALL  WRITE(A,N,N*MOUr> 

IF(ICHOICE.EO.l)  GO  TO  AC 
12  WR ITE(NCUT* 125) 

125  FORMAT !/ *1X» "ENTER  SYSTEM  B  MATRIX  (N  X  M )  :",/) 

CALL  RE  AD  !B  »ft  ,.M  ,ft  IN  ) 


G. 


CALL  WRI  TE!8  ,N,M,NOU  T> 

W« ITE!N0UT,20G) 

READ  ININ **)  IANS 
IFIIANS.EG.1)  GC  TO  12 
WR ITE I M0UT»2 1 0 ) 

210  FORMAT!/, IX, "SYSTEM  B  MATRIX  :",/) 

CALL  WRI  TE !  B  »N  »M  ,M C'J T ) 

1FUCH0ICE.LT. 9)  GO  TO  13 

11  WRITEINOUT, 127) 

127  FORMAT  </ ,  IX, "ENTER  FEEDBACK  GAIN  MATRIX  G  < M XN > I  "  . / ) 
CALL  REAQ!G,M,N*MN) 

CALL  WRI TE!G,M»N,NOUT) 

WRITEINOUT *200) 

READ  ININ  »*  )  IANS 
IFUANS.EQ.I)  GO  TO  11 
WRITE(M0UT*212) 

CALL  URITE!G,H,N,MOUT) 

1FUCHOI  CE.LT.5)  GO  TO  12 

1=  WR ITE! NOUT, 129) 

125  FORMAT!/ »1X* "ENTER  OBSERVER  GAIN  MATRIX  K  (NXM ) i  "*/) 
CALL  READ!K1,N,M.NIN) 

CALL  WR  I  TE !  K 1 ,  N  ,  M  , NOUT ) 

WRITE! NO UT *200) 

RE AO  ININ.*)  IANS 

IF  !  I  ANS. EQ  .  1 )  GO  TO  15 

URITE!H0UT*213) 

CALL  WRITE!K1,N*m,M0UT> 

GO  TC  13 

90  CONTINUE 

19  WRITE! NO UT*130) 

120  FORMAT!/,  IX,  "ENTER  SYSTEM  C  MATRIX  (M  X  N)  .'",/) 

CALL  READ!C,M*M,NIN) 

CALL  WRITE!C,M,N,NOUT) 

WRITE! NOUT, 200) 

READ  <  NIN  »* )  IANS 

IF! IANS.EQ.l)  GO  TO  19 

WRITE!MOUT*215) 

215  FORMAT!/ ,1X, "SYSTEM  C  MATRIX  J",/> 

CALL  UftITE!C,M,N,MOUT) 

12  CONTINUE 

CALL  EQUATE! A, AMAT,N  ,N) 

CALL  EQUATE ! B,BM AT, N  ,M) 

DO  5  1  =  1  ,N 
DO  5  J=1,N 

AC ! I ,J)=CMPLX(A !  I,J)  ,0.0  ) 

5  CONTINUE 

C  CALL  WRITC! AC,N,N,MOUT) 

CALL  MTR AN! A,N,N,AT) 

C  CALL  WRITE!  AT,  N,N,M0'JT) 

CO  6  1=1, V 

CO  G  J=1 ,M 

:C(I  ,J)  =  CMPLX!B! I,J) ,0.0  ) 
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o  o  r> 


6  CONTINUE 

C  CALL  WRITC(BC*N»M*MOUT) 

IF1ICH0ICE.EQ.2)  GO  TO  50 
IF!  ICH0ICE.GE.3)  GO  TO  A  1 


OBSERVABILITY  MATRIX  SECTION 


N=NSYS 

M=MC 

CALL  EQU ATE!C,CA*M»N) 

CALL  WRI  rE<CA,M*r«*,MOLT> 

CALL  WRITE!CA*M*N*NOUT> 

IT1=C 

NN-N/2 

DO  70  IT  =  1  *NN 
C 

DO  72  1=1, M 
IT1=IT1+1 
DO  72  J=1*N 
G!IT1«U>=CA(  I  *  J ) 

72  CONTINUE 
C 

CALL  MHUL!CA*A*M*N»N«CA1) 

CALL  EGUATE!CA1,CA*M*N> 

WRITE! MOUT  *220)  IT 
220  FORMAT!//. IX, 12, "TH  CA  :-,/> 

(%  CALL  URITE!CA,M,,\|, MOOT) 

7  C  CONTINUE 

WRITE!MOUT  *310) 

WRITE!M0UT,225> 

225  FORMAT!//, IX, "OBSERVABILITY  MATRIX,  Q  I",/ 
CALL  WRITE!Q.N,N»MOUTI 
CALL  EQUATE !Q»Q1,N,N) 

C 

IA  =  I  2 
IU  =  1 2 
NU=N 
C 

CALL  LSV  CF!Q,IA,N,N»UT»IU*NU,S»UK, IER) 

C 

WR I T  E  !  MO  UT  ,  A  0  0  >  IER 
C 

CO  30  1=1, N 
CO  30  J=I»N 

IFCAE?S(GCI,J)).LT.lE-05)  G! I,J>=0.0 

30  V!I, J)=Q!I*J> 

DO  31  I  =  1  *  N 
00  31  J=1,N 

IFCABS<UTlI*J)).LT*lE-05>  UT!I*J>=0.0 

31  UT 1 C I , J)  =UT ( I , J> 

WRITE! MO UT.310) 

WRITE! MO UT *230) 
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o  o  o 


a 


220  FORMAT (//*1X*"SIKGULAR  VALUES  OF  Q  J"*/) 
WR1TE(M0UT*231>  (S(X)*K=1*N) 

221  F0RMAT<1X*6E12.5,/*1X*6E12.5> 
WRITE(M0UT»23S> 

225  F0RMAT(//*1X*"UTRANSP0SE  MATRIX  :"*/> 
CALL  ««ITE<UTl,WtN».'10ur) 

CALL  MTRAN(UT1»N»N*UT) 

WRITECK0UT*236) 

226  FORMAT(//»lXt"U  MATRIX  :•%/> 

CALL  W  R 1 TE(UT  *N*N*MOLT) 

WRITE(MOUT»310) 

VIR I T  E  (MO UT  »  2  A  0  > 

240  FORMAT  <//*lX*"V  MATRIX  J"*/) 

CALL  WRITECV  *N  *N *HOU  T ) 

CALL  EQUATE <Ul»P»NtM) 

C 

GO  TO  20 
50  CONTINUE 


CONTROLLABILITY  MATRIX  SECTION 


NsNSYS 

K=MC 

CALL  EGUATE(B,3A,N*K> 

CALL  URITE(BA*N»M*MOUT> 

JT2=  0 
AN=N/2 

00  60  IT=1,NN 
C 

WRITE(M0UT*315>  IT 
315  F0RMAT(//*1X,I2,"TH  3A  :«,/) 

CALL  WRITE (BA*N*M*MOUT) 

C 

CO  61  Is 1 »N 
JT  1  =  0 

CO  61  J=1*M 
JT1=U+JT2 
P(I*UT1)=BA(  I*J) 

61  CONTINUE 
C 

CALL  MMUL<  A*BA»N*N*M«8A1) 

CALL  EGUATE<BA1«BA*N»M> 

UT2=vT2*2 
60  CONTINUE 

WRITE (MOUT  *310) 

WRITE (MOUT  *320) 

320  F0RMAT(//,1X»"C0NTR0LLA.BILITY  MATRIX*  P  : 
CALL  WR I TE ( P  *N *N  »MOU T) 

CALL  EQU  ATE ( P*P1 *N «N  > 

C 

I A  =  1  2 
IU  =  1  2 


"♦/) 
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ooo  o  ...  °  _  _  o  r>  r>  r> 


NU  =  N 

CALLING  SINGULAR  VALUE  DECOMPOSITION  ROUTINE 

CALL  LSVCF (P  *IA«N*N*UT »IU»NU  *S *WK* IER) 

WRITE<MOUT*40Q)  IER 
00  65  1=1  *N 
00  65  J=1*N 

IF<A8S(P(ItJ)).LT.lE-Q5)  P<I*J>=0.0 
£  V(IiJ):PfItJ) 

CO  66  I  =  1  *  N 
CO  66  J=1*N 

IF<A6S<UT< It J) ) .LT.1E-05  )  UT<I»U>=0.0 
6  UT1(I ,J)=UT( I»J> 

WR ITE  <M0UT«  3 10 ) 

WRITE  <MOUT*330) 

30  FQRMAT<//*1X*»SI?mGULAR  VALUES  OF  P  :"*/) 

WR ITE<M0UT«331)  <S  <K  >  »X  =  1  »N  ) 

2JL  F0RMAT<1X*6E12»5»/  ,1X,6E12.5> 

WRITE (MOOT *332) 

332  F0RMAT(//*1X*"U  TRANSPOSE  MATRIX  I"*/) 

CALL  URITE(UT1»N»N*M0UT) 

CALL  MTR4N(UT1*N*N*UT) 

WRITE(MOUT*333) 

333  FORMAT (//*1X»"U  MATRIX  J"*/) 

CALL  WRITE<UT,N*N*M0UT> 

WRITE<MOUT*310> 

WRIT£<MOUT*340 ) 

340  F0PMAT(//*1X*"V  MATRIX  J"*/> 

CALL  URITE(V»N*N.MOUT) 

CALL  EQUATE<P1*P*N«N> 

2  C  CONTINUE 

WRITE! MO UT *310) 

I0GT=3 

CALLING  LU  DECOMPOSITION  ROUTINE 

CALL  LUOATF(P*LU*N«I  A»I0GT»D1*D2  *  I P  V  T *  E fl(J  I L  »  u  A ♦  I  ER  ) 
C 

WR I TE ( MOUT * 4  00 )  IER 
CET=<01*<2**02>  > 

C 

CO  73  I=1*N 
00  73  J=1*N 
IF<J.GT.I>  GO  TO  71 
IF(J.EO.I)  GO  TG  74 
LCI *J>=LU<I *J> 

GO  TO  73 
71  CONTINUE 


73 


oooooooo  on 


Ot 


LCItJ)=0.0 
GO  TC  73 

74  continue: 

L(I»«J)  =  1.0 
72  CONTINUE 
C 

CO  80  1  =  1, N 
CO  80  d=l*N 
IF(J.LT.I)  GO  TO  61 
U(I*U)=LU(I*J) 

GO  TO  80 
81  CONTINUE 

U(I,d)=0 .0 
8  C  CONTINUE 
C 

URITE(M0UT*35G) 

350  FORMAT !//*lX ,"LU  MATRIX  :",/) 

CALL  WRITE(LU»N»N*MOUT) 

URITc!MQUT,310) 

WRITE! MO UT»352) 

352  FORMAT !/ / , IX  * "L  MATRIX  I",/) 

CALL  WRI TE!L,N,N,M0UT) 

WRITE! MO UT»310) 

WRITE (MO UT* 3 54) 

354  FORMAT ( / /, IX , "U  MATRIX  :"♦/> 

CALL  URI TE(U»N,N»M0UT) 

WRI TEC  MO  UT, 360 )  01, 02, GET 

360  FORMAT!//, IX, "01  a " , £12 . 5 ,2 X , *02  =",E12.5,2X ,"DET 
&E12.5) 

WRITE (MOUT ,362) 

362  FORMAT!//, IX, "IPVT  VECTOR  I",/) 

CALL  WRITECIPVT «N«1,MQUT  ) 

W RITE! MO UT, 364) 

364  FORMAT!//, IX, "EQUIL  VECTOR  :«,/) 

CALL  URITECEQUIL-N*! ,M0UT) 


41  CONTINUE 


STATE  TRANSFORMATION  SECTION  (X  =  TZ) 


CALCULATE  EIGENVALUES  ANO  RT  EIGENVECTORS  OF  A 

N=NSYS 
M=MC 
I  J00=2 
I A  =  1 2 
I  T=  1 2 
C 
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CALL  EIGRF  ( A  «N  , I  A*  I J OB* £  VAL  *  T*  I  T  *  JK  1 .1 ER > 


CO  85  1  =  1  «N 
£R<I)=R£AL!EVAL< I) ) 

IF<ABS<ER< I) ).LT.lE-06>  ER  <  I  )=0 . 0 
£1  Cl  )=AIMAG<£VAL  <1 ) ) 

IF<ABS<EI<D  ).LT. IE- j6>  El  CD  =0.0 
CVAL!I)  =  CMPLX<ER!I)»*I!  I  )) 

EVMAG<D  =  <ER  <1)**2  ♦  El  <  I )  >  *  *  .5 

CONTINUE 

W»ITE<MOUT*310) 

WRITE<MOUT»400)  IER 
FORMAT  </ ,1  X, "IER  =",13) 

UR  IT  E  <M0UT ♦ A  02 ) 

FORM  AT <//, IX, "EIGENVALUES  J "  ,// , 6X , "REAL ",8X , « IH AG" , /) 

CO  83  1=1, N 

URITE<MOUT»404)  £VAL<I> 

CONTINUE 

FORMAT  <1X,2E12»5> 

URIT£!MOUT,406) 

FORMAT!//, IX, "EIGENVECTORS  <COMPLEX  TR ANSF OR  MAT  1 CN  MATRIX)' 
S"  I" »//,5X,3<"REAL",8X,"  IMAG",8X),/ ) 

CALL  CLEANC<T,N,N> 

CALL  URITC(T,N,N,MOUT) 

URIT£<M0UT,411)  WK1<1) 

F0RMAT</*1X, "PERFORMANCE  INDEX  =",E12.5> 

FORM  REAL  TRANSFORMATION  MATRIX 

CALL  FCRM<T*N,N,TREAL> 

CALL  EQUATE! TREAL,TLU,NtN) 

CALL  MTRAN<TREAL*N,N,TTRAN) 

*«R ITE  <MOUT, 4  08 ) 

FORMAT!//, IX, "TRANSFORMATION  MATRIX  T  (READ  :",/) 

CALL  URIT£(TREAL,N,N,MOUT) 

URITE<MOUT,310) 

CALL  MMUL<TTRAN,TR£AL,N,N,N,TTT) 

CALL  CLEANK<  TTT,N,N) 

WRITE!. MOUT, 407) 

FORMAT!//, IX, "TTR AN *T  MATRIX  1",/) 

CALL  URITE(TTT*N,N,MOUT> 

WRITE! MO UT, 310) 

FORM  9X9  SU9 -TRANS FORM AT  ION  MATRIX 
11=0 

CO  35  1=1, N 
IPCI.LE.31  GO  TO  35 
11=11*1 
J1  =  0 


nooo  r>  o  ono 


CO  3b  J=  1  «N 

IF ( J  «LE  • 5 )  GO  TO  36 

t 1 — J 1+1  • 

T9X9(Ilt  Jl)=TREAi.(  I»  J) 

3b  CONTINUE 
35  CONTINUE 
N9  =  1 1 

CALL  URITE(T9X9«NS,N9,H0UT) 

CALL  E3UATE<T9X9«T9,U9tN3) 

IA-12 

IDGT=0 

C 

CALL  LINV2F(T9,N9»IA*T9I,I0GT»yK3» IER) 

C 

URITE(M0UT*419) 

A  19  FORMAT (//tlX*"T (9X9)  INVERSE  I  "»  /) 

CALL  WRITE<T9I»N9tN9*M0UT> 

CALL  MMUL(T9I,T9X9»N9, N9»N9,TTT9) 

CALL  URI TE<TTT9«N9*NS*M0UT> 

CALL  MTR AN(T9X9»N9*N9*TT9X9) 

CALL  MMUL(TT9X9*T9X9«N9*N9.N9,TTT9) 

URITE(MOUT»310) 

URITE<MOUT*409) 

409  FORMAT <//«lX»"TTRAN  19X9)  *  TC9X9)  MATRIX  I",/) 

CALL  CLE ANR1TTT9»N9«N9) 

CALL  URITE(TTT9*!\i9*N9»M0UT> 

CALCULATE  INVERSE  OF  THE  REAL  TRANSFORMATION  MATRIX  T 

I  A  =  1 2 
IDGT=0 

CALL  LINV2F (TLU»N» IA*TINVR» I DGT  ♦ 'JK3  » IER ) 
yRITECMOUTtAOO)  IER 

CALCULATE  THE  INVERSE  OF  THE  COMPLEX  TRANSFORMATION 
MATRIX  T  (A  COLUMN  AT  A  TIME) 

M=N 

IJOB=l 
C 

CALL  LE02C(T»N*12»TIVVC»Mtl2»IJ02»WA»WK2»IEP ) 

WRI TECMOUT »4  00 )  IER 
C 

MSI 

I  J0B  =  2 
CO  B 7  Js 1 «N 
DO  88  I  =  X t N 
IF(I.EQ.J)  GO  TO  89 
B2(I)sCMPLX<0,0»0.0) 

GO  TO  88 
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65  CONTINUE 

6Z!I>=CMPLX(1.0»0.C> 

Pi  CONTINUE 

CALL  LEQ2CiT»N»12,BZ,M»12,IJOB,<.A*wK2*IEK> 
WRITE(MOUT,400)  IER 
DO  90  1=1, N 
TINVC!I»J)=BZ(I> 

9  C  CONTINUE 
87  CONTINUE 

yAITC(M0UT,412) 

♦12  FORMAT!//, IX, "T  INVERSE  MATRIX  ( COMPLEX)  :",/> 

CALL  CLEANC (TINVC,N,N) 

CALL  WRI TC( riNVC»N»N»MOUT) 
yRITE(MOUT,310) 

WRIT£!M0UT»414) 

♦14  FORMAT!//, 1X»*T  INVERSE  MATRIX  (REAL)  :"♦/> 

CALL  CLEANR (T1NVR,N,N) 

CALL  WRITE(TINVR,N,N,MOUT> 

IF ( JP  ASS  .£Q » 1 )  GO  TO  99 

calculate  left  eigenvectors  of  a 

*****  RIGHT  NOW  THIS  IS  PASSED  OVER,  BUT  IT  WORKS  ***** 

WR1TE(MOUT,310) 

I JOB  =2 
IT  =  12 

CALL  EIGRF(AT,N,12,I00B,EVALT,TT  ,IT,«K1, IER) 

C 

DO  93  1=1, N 
ETRCI)=REAL(EVALT<  I)  ) 

IFCABSlETRd  )).LT.1E-0S)  ETR!I>=0.0 
ETI(I)=AIMAG<EVALT(I  )) 

IFtABSCETI  (I  )).LT.lE-06)  ETUI)  =  0.0 
EVALTCI)=CMPLX(ETR(I),ETI(I)  ) 

9  2  CONTINUE 
C 

WR ITE<  MOUT, 4  00 )  IER 
WRI TE (MO UT, 402) 

DO  91  1=1, N 

VIRITE(M0UT,4C4)  EVALT(I) 

91  CONTINUE 

WRITE<MQUT,406) 

CALL  CLEANC  <  TT,N  »N ) 

CALL  WRI TC!TT,N,N,MOUT) 

URITE(M0UT*4ll>  WK1(1) 

99  CONTINUE 
C 

■R I Tt (MOUT ,3 1 0 ) 

CALL  MMULC(TINVC,T «nl  ,N,N  «T1  D ) 

WRITE! MOUT, 430) 

*»J9  FORMAT!//,  IX,  "TINV*T  MATRIX  (COMPLEX)  I",/) 
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CALL  CLE  ANC ( T ID ♦ N  »N ) 

CALL  WHITC(TID*N*N*MOUT) 

CALL  MMULC (AC*T«N*N*.\*T1> 

CALL  CLEANC( T1*N*N> 

C  CALL  WR I TC ( T 1 »N  *N*MOUT ) 

CALL  MMULC(TINVC.*UtN*N*Nf2> 

CALL  CLE  ANC  ( Z*N*N) 

C 

WRITE(M0UT»31Q) 

WRITE(MOUT.420> 

420  F0RMAT(//*1X*"DI  AGONALIZED  A  MATRIX  (COMPLEX)  :«.//*?X* 
&2<"R£AL" ,8X**IMAG".8X>*/) 

CALL  URITC(Z*N»N*MOUT) 

C 

K  =  MC 

CALL  MMULC (TINVC*BC*N«N»M«TINVB) 

WRITE(MOUT*310> 

URIT£(M0UT»422) 

422  FORMAT  (//*IX»"TINV*f3  MATRIX  (COMPLEX)  .:"♦/> 

CALL  WRITC(TINVB»N*M»MOUT) 

C 

WRITE(MOUT*310) 

CALL  MMUL(TINVR,TREAL*N* N*N*TRID) 

WRITE (MOUT. 432) 

432  F0RMAT(//,1X,"TINVR*TREAL  MATRIX  :“*/> 

CALL  CLE ANR ( TR I D  *N  *N ) 

CALL  WRITE(TRID,N*N*MOUT) 

WRITE<M0UT,31C> 

CALL  MMUL(A*TREAL»N*N»N*ATR) 

CALL  CLEAN* (ATR*N*N) 

C  CALL  URITE(ATR,N,N*MCUT) 

CALL  MMUL(TINVR«ATR*N*N*N»ZREAL) 

CALL  CLEANS(ZREAL*N*N) 

WRITE(M0UT,424) 

424  FORMAT (/ / *1X*"BLGCK  DIAGONALIZED  A  MATRIX  (REAL)  :«*/> 
CALL  URITE(ZREAL*N,N»MOUT) 

CALL  MMUL(TINVR*B»N»-\*MC»TINVBR> 

WRITE! MO UT *426) 

42b  F0RMAT(//*1X»"TINV*8  MATRIX  (REAL)  :**/> 

CALL  URITE(T  INVBR *N ♦ MC ,MOUT  ) 


CONTROLLER  CLOSED  LOOP  CALCULATIONS 


IF(ICHGICE.LT.4>  GO  TO  52 

FORM  9X9  SUB -SYS TEM  COMPONENT  MATRICES  FROM 
FULL  SYSTEM  MATRICES 


11=0 

CO  23  1=1, N 
Ir  ( I  »L£  •  3 )  GO  TO  23 
11=11*1 
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CO  24  J=1*M 
H9<IltJ)=TIMVBS( I,  J) 

24  CONTINUE, 

22  continue: 

c 

CO  26  I  =  lt.M 
J1=0 

CO  27  J=1*N 
1F(J.LE.3>  GO  TO  27 
wl=Jl*l 

G9U  *U1> -G(  I«J) 

27  CONTINUE 
26  CONTINUE 

C 

11  =  0 

DO  28  1  =  1, N 
IFCI.LE.3)  GC  TO  28 
11=I1*1 
01=0 

DO  29  J=ltN 
IFTJ.LE.3)  GO  TO  29 
J1=J1+1 

A9TI 1,J1 )=2REAL(I,J) 

25  CONTINUE 

28  CONTINUE 

PERFORM  FUL  12X12  CLOSED  LOOP  SYSTEM  ANALYSIS 

CALL  MMUL<TINVBfifG,N,MfN,9G) 

URITETMOUT  ,310) 

CALL  WRITECBG,N,N,r10UT) 

CALL  MAOO<ZREAL,BG,N,N,l  .OtABG) 
fcRITETMOUTf 310) 

WRITEOIOUT,  440) 

440  FORMATC//* lXt "DIAGC A)  ♦  EG  MATRIX  1"*/) 

CALL  CLEANR ( ABGf N»N) 

CALL  WRITE(A8GfN*N*M0UT) 

C 

IJOB=0 
1 ABG=12 
I T  =  1 2 
C 

CALL  EIGRF(ABGfN*IABGfIJ03fEVALTfTTfIT»WK4fIER) 
C 

CO  54  1=1, M 
ETR(I)=REAHEVALT(I)  ) 

IFTABSCETRC 1) ) .LT.1E-06)  ETR<I)=0.0 
ETI ( I)=AIMAG (EVALTt I )) 

1F<ABS(ETI(1) ) .LT.1E-06  >  ETI(I)  =  0.0 
EVALT(I)=CMPLX<ETR(I ),ETJ(l) ) 

54  CONTINUE 

WRITE(MOUT»310) 
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WRITE!MOUT,4GO>  IER 
WRITE!M0UT,4Q2> 

CO  S5  I- 1*N 

yRITE!M0UT,4C4>  £V  AL  T ( I > 

55  CONTINUE 

PLUGGING  CONTROL  BACK  INTO  ORIGINAL  NON-DIMENSIGNALI ZED 
SYSTEM 

WRITE ! MO  U  T  *  3 1 C ) 

CALL  ESUATE!  AMAT,A».\I,N) 

CALL  E3U ATE (RMAT  »B,N  *M) 

CALL  WRtTE!A,N,N»MOUT) 

CALL  WRITE!8»N*M,M0UT) 

WR IT  E (MOUT*  3 1 0 ) 

CALL  MMUL! 9  »  G  ?  N  »  M  »  N  ,  BG  ) 

CALL  URITE<BC-,N,N,M0UT> 

CALL  MMUL!3G,TINVR,N,N,N,BGTI> 

CALL  WRITE !8GTI,N,N»M0UT  ) 

CALL  MADD! A«  fiGTI «N,N,1.0  * ABG  ) 

WRITt!MOUT,310> 

WRITE! MOUT *442) 

442  FORMAT!//, IX, "A  ♦  8*G*TIAV  MATRIX  :-,/> 

CALL  CLEANR! AOG,N,N> 

CALL  WRITE!ABG*N,N«MCUT) 

I JOO-O 
1ABG=12 
IT  =  12 

CALL  EIGRF!ABG,N,IA3G,IJCB,EVALT,TT,IT,WK4»IER> 

CO  5b  1=1, N 

ETRt  I)=REAL!EVALTm  ) 

IF!A3S!ETR!I )).LT.l£-06)  ETR!I>=0,0 
ETI!  I)=AIMAG!EVALT!D) 

IF!AES!ETl!I)).LT.lE-06)  ETI !I>=C.O 
EVALT(1>=CMPLX!ETR!II*ETI!I> > 

56  CONTINUE 
WRITE!MOUT,310) 

WR1TE!MQUT,400)  IER 
WRITE!MOUT,402) 

CO  57  1=1, H 

WRITE!MOUT,404)  EV  AL  T 1 1 ) 

57  CONTINUE 

PERFORM  9X9  CLOSED  LOOP  SYSTEM  ANALYSIS 

WRITE!MOUT*310) 

CALL  WRITE! A9,N9,N9, MOUT  ) 

CALL  WRITE!B9,N9,M,MGUT) 

CALL  WRI TE!G9,M,n9,M0UT) 
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CALL  MMUL(89,G9,n9,N9,N9,8G9) 

CALL  WRI TE(BG9,N5,N9,M0U f) 

CALL  MA00(A9,BG9,N9»A9*1  ,0»ABG9) 

WRITE! MO UT ,460) 

FORMAT!//, 1X»"DIAG{A9)  ♦  BG9  MATRIX 
CALL  CLEANH <ABG9,N9,N9) 

CALL  MR I TE ( ABG9 , N9  »N9,M0LT) 

IJ09=0 
I  ABG=12 
IT=12 

CALL  EIGKF(A3G9,N9»IABG, IJ03,EVALT,TT,IT ,WK4,IER  ) 

WRITE(MOUT,400>  IER 
WRITE!MOUT,402> 

00  62  1= 1,N9 

WRITE<MOUT,404>  £VALT<I) 

CONTINUE 

CALL  WRITE(T9X9»N9»N9»M0LT) 

CALL  WRITE(T9I,.\S9,N9,M0UT) 

CALL  MMUL!ABG9*T9X9,N9*N9,N9,A91> 

CALL  MMUL<T9  I,A91«N9»N9*N’9«A9CL) 

WR ITE (MO  UT  «  4  62 ) 

FORMAT(//,lX, "ORIGINAL  SYSTEM  CLOSED  LOOP  A9  MATRIX!-,/) 
CALL  CLEANR! A9CL,N9,N9) 

CALL  WRITE!A9CL,N9*N9»M0UT> 

IJ06-0 
IABG=12 
IT  =  12 

CALL  EIGRFCA9CL,N9»IABG,IJ0S,EVALT,TT,IT,WK4,IER  > 

WRITE<MOUT,400)  IER 
WRITE(MOUT*402) 

DO  63  1=1, N9 

WR  ITE!MOUT*404)  EVALT(I) 

CONTINUE 


C!l,l>=-. 001329 
C(l,2)=* 03043 
C( 1 ,3 )=• 603  4 

C<  2 *  1 ) =• 7068  3 
C(2,2)=. 08840 
C(2,3)=. 68922 
WRITE (MOOT ,310) 

CALL  WRITE(C,M,N,MOUT) 


CALL  MMUL(C,TREAL,M,fi,N,CTREAL> 
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yRIT£(M0UT.444) 

444  F0RMAT(//.1X."C*TREAL  MATRIX  I"./) 

CALL  URITEtCTPEAL.M.N.MOL’T) 

IF<ICH0ICE.LT.5>  6C  TO  52 
C* **********************************  *********** 

C  OBSERVER  CLOSED  LOOP  CALCULATIONS 

O ************************************  ********* 

VR I TE(M0UT«  3 1 0 ) 

CALL  MTR AN<CTR£ALtM.NtCTT) 

CALL  MTRANtKl.iJ.M.KlT) 

C  CALL  WRITEtCTT.N.M.MGUT) 

C  CALL  WRITE(KIT.M.N.MOUT) 

CALL  MMUL ( CT  T  »K  1  T  »N»P»N*i<CT  > 

C  CALL  MMULCK 1 1 CTR EAL  tN  til  t  Nt KC  T) 

CAUL  CLEANS  (KCT.NtN) 

CALL  URITE(KCTtN.N.MOUT) 

CALL  MADO(ZREAL.KCT. N.N. 1.0.AE2) 

WRITE<M0UT.448> 

448  FORMAT(//tlX  ."OBSERVER  (LAMDA  -  CT*KT)  MATRIX  i",/) 

CALL  URITE(AE2.N.NtM0UT) 

ALL  EQUATE (AE2.AE2A.NtN) 

3 

uoa=o 

IA£= 12 
T=  12 
C 

CALL  EI(5RF(  AE2.NtIAE.lJC3tEVALT.TT.ITfJK4t  IFR) 

C 

WR ITE ( MOUT*  4  02) 

CO  25  I  =  1.N 

WRITE<M0UT.404)  EVALT(I) 

25  CONTINUE 

CALL  EQUATE* AE2A.AE2. N.N) 

C  CALL  yRITE(AE2.N.NtM0UT) 

CALL  MMULCAE2.TINVR.N.N.Nt AE1) 

CALL  MMUL(TREAL.AEl.N.N.NtAE) 

WRITECMOUTt 310) 

WRITE(MQUT*450) 

450  F0RMAT(//,1X."CL0SED  LOOP  ERROR  STATE  COEFFICIENT  MATRIX 

&-  :"./) 

CALL  CLEANR* AE.N.N) 

CALL  WRITE(AE.N.NtMCUT) 

C 

IJOS=Q 
IAE=12 
I  T=  1 2 
C 

CALL  FIGRFtAE.N.IAE. JJ0B.EVALTtTTtIT.hK4.IER) 

C 

CO  21  Is 1 ♦ N 

ETR( I)=REAL( EVALT( I ) ) 

1F(A2S(ETR ( I ) > .LT.1E-06)  ETR(I)=0.0 


oonooo  o  o  o  r»  o o  r>  o  o  n 


ETI  < I)=AIMA6(EVAlT(I>) 

IF(ABS(£TI (I > > .LT.1E-06 >  ETI  (I)  =  Q.O 
EVALT(I>=CMPLX(ETR(n,ETI(I>  > 

21  CONTINUE 
4RITE(MOuT,400>  IER 
*RITE(MOUT,402> 

00  22  I=1*N 

UR  IT  E (MOUT, 4  G4>  EVALT(I) 

22  CONTINUE 


52  CONTINUE 
STOP 
END 

*«********,«*»*******»*******,****.*,**** 

SUBROUTINE  READF (MAT ,N » M  ,L IN > 

A**************#*********.**********.**** 

THIS  SUBROUTINE  READS  THE  REAL  STATE  MATRICES 
FROM  DEVICE  NUMBER  LIN 

REAL  MAT (12 *  12 ) 

DO  1  LS1 *M*S 
K=L*5 

1F(M-L.LT.5)  K=M 
CO  2  1  =  1, N 

REAO(LIN.IOO)  (MAT<I,J>,U=L,K> 

2  CONTINUE 

ICO  FORMAT (SE12.5) 

REAO(LINtllO)  XXX, YYY 
110  F3RMAT(A1*/*A1) 

1  CONTINUE 

RETURN 

END 


SUBROUTINE  CLE AMR ( M A T* N , M> 


THIS  SUBROUTINE  ZEROES  OLT  SMALL  ELEMENTS  OF 
A  REAL  MATRIX  MAT ( N, M ) 

REAL  MAT (12,12) 

C 

DO  l  1  =  1  ,N 
DO  1  J=1  ,M 

IF  (ABS (MAT (I *J) ) .LT. IE-0  5)  MAT(I,J)  =  0.0 
1  CONTINUE 
C 

RETURN 

CNU 
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C* ************************  ******* ********* 

SUSR  GUT  I  N£  CLEANS(MAT,N,M> 

C* ****************************  ************ 

v 

C 

REAL  M'AT  (12*12) 

C 

CO  1  1=1 «N 
00  1  J=1 ,M 

IF (AES (MAT (I  ,J)  )  .LT. IE-02)  HATH  *J)  =  0.0 
1  CONTINUE 
C 

RETURN 

END 

C****************************  *************** 

SUBROUTINE  F ORM < A ,N, M.3 > 

THIS  SUBROUTINE  FORMS  A  REAL  TRANSFORMATION  MATRIX 
FROM  A  COMPLEX  TRANSFORMATION  MATRIX 

REAL  AR(12,12)«AI(12«12) ,8(12,12) 

COMPLEX  A( 12 , 12 ) 

C 

DO  1  1=1, N 
DO  1  J=1  ,M 
AR(I,J>  =  REAL(A(I,J>> 

AKI  ,J)  =  AIMAG(AU,J>) 

1  CONTINUE 
01=1 

CO  3  J=1  »N 
DO  2  1  =  1  ,N 

IFIABS(AKItJl))  .LT.1E-10)  GO  TO  2 
GO  TO  5 

2  CONTINUE 
CO  8  1  =  1, N 

8  ( 1 , 01  )  =  AR  ( I  *  J1 ) 

8  CONTINUE 

IF(Jl.GE.N)  GO  TO  20 
Jl=Ul*l 
C-0  TO  3 
5  CONTINUE 

CO  10  1=1, N 
B(l«ol)=AR(I«01) 

6(1,01*1 )=AI (1,01) 

1C  CONTINUE 

lF(Jl.GE.N)  GO  TO  20 
ulsjl+2 

3  CONTINUE 

2  C  CONTINUE 

C 

RETURN 

END 


oooooo  nnoooo  n  r>  ooooo 


SUBR  OUTINE  C LE AN C ( M A  T ,  N •  M ) 


THIS  SUBROUTINE  ZEROES  OUT  SMALL  ELEMENTS  OF 
A  COMPLEX  MATRIX  MAT  (N  «M  ) 

REAL  MAT R (12*12) .MATIC12.12) 

COMPLEX  MAT(12*12) 

CO  1  1=1. N 
DO  1  J=1*M 

MATR  (I,J  )=REAL(HAT(I  ,J)  > 

IF (A8S(MATR ( I  *  J ) )  •  LT  »1E “ C6 )  M A TR ( I ,J  )  =  0  .  C 
MAT  I (  I«J)=AIMAG(MAT ( I  •  J  ) ) 

IF(A6S(MATI(I«J) ).LT.lE-06>  M AT  I (I  * J  )= 0 . 0 
MAT( I* J)=CMPLX(MATR( I* J >  *M A T I ( I . d) ) 

1  CONTINUE 

RETURN 

END 


SU5RCUTINE  EQUATE(A*3*N1 ,M1> 
***************************************** 

THIS  SUBROUTINE  EQUATES  REAL  MATRIX  B  TO  REAL 
MATRIX  B  (A  IS  INPUT*  3  IS  OUTPUT) 

REAL  A(12*12). 6(12*12) 

00  1  1=1 *N1 
CO  1  J=1,M1 
3(I*d)=A(I«d) 

1  CONTINUE 
RETURN 
END 

************************* ************ ******* 

SUBROUTINE  E3UATC ( A* c»Ml *M1 > 


THIS  SUBROUTINE  EQUATES  COMPLEX  MATRIX  3  TO  COMPLEX 
MATRIX  A  (A  IS  INPUT,  6  IS  OUTPUT) 

COMPLEX  4(12,12) *B(12,12) 

DO  1  1  =  1  »M1 
DO  1  J=1*M1 
t(I*d)=A(I*d) 

1  CONTINUE 
RETURN 
END 

C* **********************************  ****** 

SUBROUTINE  M  MUL ( X , Y , N 1 ♦ N  2  *N  3  *  Z ) 

C* **********************************  ****** 


no  o  oooooo  o  ooooo 

»  *  m  i\>  m  »  »  cw  fo  H 


THIS  SUBROUTINE  MULTIPLIES  TWO  REAL  MATRICES  J 


XCN1XN2)  X  Y (N2XN3 )  =  2(N1XN3> 

REAL  X<12,12)*Y<12»12>*Z<12*12) 
CO  3  J=1*N3 
CO  2  1=1 *N1 
S=0. 

CO  1  K=1 *N2 
S=S*X<I»K)*Y (K  *U ) 

2<I*U)=S 

CONTINUE 

RETURN 

ENO 


SUBROUTINE  MMULC <X ♦ Y *N1 *N2 *N3» 2 > 


THIS  SUBROUTINE  MULTIPLIES  TWO  COMPLEX  MATRICES 

X (N1XN2)  X  Y (N2XN3  )  =  2(N1XN3> 

COMPLEX  XC 12 *12 >*Y< 12,12  >*2(12.12) *S 
DO  3  J=1 «N3 
CO  2  I=l.Nl 
$=CMPLX(0. 0*0.0) 

CO  1  K=1 »N2 
S=S«-X< I*  K)  *  Y ( K  *  J ) 

2CI*U)=S 

CONTINUE 

RETURN 

END 


SUBROUTINE  MTRAN ( A*NR»NC  *B> 


REAL  A(12*12)«8( 12*12) 

C 

00  1  1=1 ,NR 
CO  1  J=1»NC 
1  E ( J *  I ) =A ( I « J  ) 

C 

RETURN 

ENO 

C* ***...**••*•**»+.*«**•. ********* 
SUBROUTINE  M  T  R  AN  C ( A  *  NR  *  N  C  *  0  ) 
C* ******  ******************  ******** 

c 

COMPLEX  A(12tl2) ,3(12,12) 

C 

00.1  1=1 »NR 
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CO  1  J=1»NC 
1  E(J.I>=A(I,J) 

C 

RETURN 

END 

C#  »**«*******»*•**»••«*•**•**•  ******** 
SUBROUTINE  MADO < A, 3, N. M. CONS T. C> 


THIS  SUBROUTINE  AODS  TWO  REAL  MATRICES  *. 

A(NXM)  ♦  CONST *B (N  XH )  =  C(NXM) 

REAL  A(12»12) .6(12.12) »C (12.12) 

CO  1  1=1 «N 
CO  1  J=1 tM 

C(  I.O)=A(I, J> ♦CONST* o< I • J> 

1  CONTINUE 
END 


SUBROUTINE  READ ( Y.N. K. J I N) 


THIS  SUBROUTINE  IS  USED  TO  READ  IN  A  REAL  MATRIX 
Y(NXM)  FROM  DEVICE  NUMBER  JIN 

REAL  Y<12«12> 

NM=N*H 

1C  CONTINUE 

READ (JIN  »*  )  I tJ* VALUE 
IF(I.EQ.O)  RETURN 
Y(I.J)=V ALUE 
60  TO  13 
END 


SUBROUTINE  WRITE<MAT  *N*M  »JOUT) 


THIS  SUBROUTINE  IS  USED  TO  WRITE  A  REAL 
MATRIX  MAT(NXM)  TO  DEVICE  NUMBER  JOUT 

REAL  MAT (12*12) 

CO  1  L=1.M*6 
X=L*5 

IF(M-L*LT.5>  K=M 
CO  2  1=1. N 

WRITE! JOUT. 100)  (MAT ( I  * J >,J  =  L»K  ) 
i  CONTINUE 

100  FORMAT (6E1 2 • 5 ) 

WRITE(JOUT.UO) 

110  FORMAT (//) 

1  CONTINUE 
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C  ytITECdOUT.no) 

RETURN 

END 

C *«*.*«*•*•*** ••*****•. **»•♦. ******* 
SUBROUTINE  WRITC(HAT .N.M.dOUT) 


THIS  SUBROUTINE  IS  USED  TO  WRITE  A  COMPLEX 
MATRIX  MAT(NXM)  TO  DEVICE  NUMBER  dOUT) 

REAL  MATR(12*12)»MATI(12»12) 

COMPLEX  PAT (12.12) 

C 

00  1  L-l.M.d 
K~L*2 

IF(M-L*LT.2)  K-M 
CO  2  1  =  1  .N 

URITECdOUT.lOO)  (MAT < I .d > .JaL.K ) 

2  CONTINUE 
100  FORMAT (6E12» 5) 

.RITE(dOUr.llO) 

110  FORMAT (//) 

1  CONTINUE 
C  WRITE (dOUT* 1 10 ) 

RETURN 

END 
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